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Abstract 


A  new  time-domain  method  for  the  analysis  of  planar  guided  wave 
propagation  and  scattering  is  developed  in  which  an  analytical  process  is  incoporated 
along  one  of  the  spatial  dimensions.  The  method  makes  use  of  the  Method  of  Lines 
which  has  been  applied  to  microwave  problems  in  the  time-harmonic  cases. 

The  procedure  and  the  formulation  of  the  proposed  method  to  calculate  the 
time  history  of  a  pulse  scattered  in  a  partially  filled  waveguide  and  a  finned  waveguide 
are  described  in  detail.  The  cutoff  characteristics  of  the  structures  can  be  derived  via 
the  Fourier  transform  of  the  time-domain  pulse  scattering  data.  The  results  are 
compared  with  other  available  data. 

The  method  is  extended  to  characterize  three-dimensional  structures.  The 
procedure  and  the  formulation  of  the  three-dimensional  analysis  are  explained.  The 
propagation  and  scattering  data  of  a  pulse  in  a  planar  transmission  line  with  its 
discontinuities  can  be  obtained  by  using  the  three-dimensional  analysis.  From  the 
time-domain  data,  the  frequency-domain  characteristics  for  a  wide  range  of 
frequencies  can  be  extracted  by  the  Fourier  transform.  The  important  design 
parameters  such  as  the  characteristic  impedance  and  the  propagation  constant  of  a 
uniform  microstrip  line  and  the  scattering  parameters  of  its  discontinuities  (step-in- 
width,  open-end  and  gap  discontinuities)  are  presented  and  con.;  art  r  with  available 
published  data. 
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Chapter  1 
Introduction 


1.1.  Background  and  Goals  of  This  Dissertation 

A  discontinuity  in  a  planar  transmission  line  circuit  is  caused  by  an  abrupt 
change  in  the  geometry  of  the  strip  conductor,  as  shown  in  Fig.  1.1.  Therefore, 
electric  and  magnedc  field  distributions  are  modified  near  this  discontinuity.  At  low 
frequencies,  the  discontinuity  dimensions  are  much  smaller  than  the  wavelength.  In 
this  case,  the  altered  electric  and  magnetic  field  distribution  can  be  approximated  by 
lumped  element  equivalent  circuits.  However,  as  the  operating  frequency  increases, 
accurate  frequency-dependent  scattering  matrix  representations  associated  with 
discontinuities  are  necessary  for  a  accurate  design  of  a  microwave  and  millimeter- 
wave  circuit.  Especially,  the  present  cut-and-try  cycles  in  the  Microwave  Integrated 
Circuit  (MIC)  /  Monolithic  Millimeter-Wave  Integrated  Circuit  (MMIC)  will  be  gready 
reduced  if  the  accurate  frequency-dependent  characteristics  of  discontinuity  can  be 
obtained. 

In  the  early  stages  of  the  study  of  discontinuities,  analyses  were  mostly  done 
by  quasi-static  methods  [1-10],  In  all  these  methods,  the  following  assumptions  are 
implied:  (i)  the  size  of  the  discondnuity  is  small  compared  with  the  wavelength  so  that 
the  phase  variation  across  the  discontinuity  can  be  neglected,  (ii)  the  current  on  the 
strip  has  no  divergence  and  (iii)  the  strip  conductor  is  infinitely  thin. 
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Fig.  1.1(a)  The  layout  of  a  microwave  amplifier  using  a  gallium  arsennide  MESFET 
device,  showing  discontinuities  in  the  microstrip  lines. 


Eig-  1.1(b)  Various  micro  strip  discontinuities  appearing  in  microwave  integrated 
circuits,  (i)  open-end,  (ii)  gap,  (iii)  step-in-width,  (iv)  T-junction, 
(v)  cross- junction  and  (vi)  bend 
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Because  of  these  assumptions,  quasi-static  analysis  is  valid  only  up  to  a  few 
gigahertz.  Therefore,  fullwave  analysis  is  required  to  determine  the  frequency- 
dependence  of  various  parameters  needed  for  a  complete  characterization  of  the 
discontinuity.  The  first  fullwave  analysis  was  done  by  using  a  planar  waveguide 
model  [11,12],  where  the  frequency  dependent  scattering  parameters  can  be  obtained 
by  matching  the  field  in  each  waveguide  at  the  discontinuity  interface  using  the  planar 
waveguide  model.  The  Spectral  Domain  Method  /  Transverse  Resonance  Technique 
can  also  be  used  to  obtain  the  frequency-dependent  scattering  parameters  for  various 
discontinuities  [13-16],  More  recently,  the  moment  method  has  been  used  by  several 
investigators  to  characterize  discontinuities,  where  radiation  and  surface  wave  effects 
can  be  included  [17,18], 

All  of  the  above  mentioned  approaches  have  been  performed  in  the  frequency 
domain;  that  is,  the  data  for  the  whole  frequency  range  are  calculated  one  frequency  at 
a  time.  This  is  an  expensive  aoproach  when  the  results  over  a  wide  frequency  range 
are  sought.  Since  a  pulse  response  contains  all  the  information  of  a  system  for  the 
whole  frequency  range,  it  is  natural  to  use  a  pulse  to  excite  a  planar  circuit  in  a  time- 
domain  approach.  From  the  time-domain  pulse  response,  the  frequency-domain 
characteristics  of  the  structure  can  be  extracted  via  the  Fourier  transform  Therefore, 
the  time-domain  analysis  of  microwave  planar  transmission  structures  provides  an 
alternative  to  time-harmonic  approaches,  and  it  is  also  useful  for  studying  the 
behavior  of  pulsed  signal  in  such  structures  as  high  speed  digital  circuits.  The 
knowledge  of  time-domain  propagation  and  scattering  can  be  used  for  circuit 
diagnosis  by  means  of  a  fast  pulse. 
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A  typical  time-domain  analysis  requires  a  discretization  of  the  three- 
dimensional  space  into  the  three-dimersional  mesh,  as  shown  in  Fig.  1.2.  Maxwell's 


Fig.  1.2.  Examples  of  the  discretization  scheme  in  (a)  TLM  method  in  two-dimensional 
plane,  and  (b)  FDTD  in  three-dimensional  space. 

equations  can  be  discretized  at  these  points,  as  in  the  case  of  Finite-Difference  Time- 
Domain  method  (FDTD)  [19],  or  the  wave  phenomena  can  be  modeled  by  networks, 
as  in  the  case  of  the  Transmission  Line  Matrix  (TLM)  method  [20]  and  the  Bergeron's 
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method  [21]  These  techniques  have  been  applied  by  several  investigators  to 
microwave  planar  transmission  lines  in  order  to  study  the  time-domain  pulse 
propagation  and  scattering  phenomena,  and  the  frequency-domain  data  have  been 
derived  from  the  time-domain  information.  In  [22]  and  [23],  Bergeron’s  method  was 
used  to  obtain  the  qualitative  results  that  graphically  illustrate  the  pulse  propagation 
and  coupling  in  the  microstrip-bend  and  the  crossing.  In  [24],  the  TL\1  method  was 
used  to  obtain  some  characteristics  of  fin  lines.  In  [25]  and  [26],  the  FDTD  method 
has  been  used  to  obtain  the  behavior  of  a  pulsed  signal  along  a  uniform  microstrip  line 
with  various  discontinuities,  from  which  the  frequency-dependent  design  information 
for  the  uniform  microstrip  line  and  its  discontinuities  was  extracted.  Usually,  large 
amount  of  computer  storage  and  long  computation  time  are 


Fig.  1.3.  Several  planar  transmission  lines  used  in  MICs.  (a)  microstrip,  (b)  slotline, 
(c)  coplanar  waveguide,  and  (d)  coplanar  strip. 
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required  for  these  aforementioned  time-domain  analyses.  An  additional  problem  of 
these  methods  is  the  difficulty  in  handling  open  boundaries. 

In  this  study,  a  new  time-domain  method  is  proposed,  that  originates  from  the 
fact  that  most  of  tne  discontinuities  in  the  planar  transmission  structure  are  loc  ted  on 
the  substrate  surface,  and  the  space  above  and  below  this  surface  is  uniform  and 
homogeneous  as  shown  in  Fig.  1.3.  One  wants  to  make  use  of  this  structural  feature 
and  wishes  to  solve  tne  problem  by  discretizing  the  structure  in  only  the  two- 
dimensional  surface  on  the  substrate  instead  of  using  the  conventional  three- 
dimensional  discretization.  This  is  possible  if  the  wave  scattering  phenomena  in  the 
direction  normal  to  the  substrate  surface  can  be  derived  analytically  By  using  an 
analytical  solution  in  one  direction,  the  dimension  of  the  problem  is  effectively 
reduced  by  one  so  that  computer  storage  and  computation  time  may  be  saved.  The 
present  method  actually  incoporates  this  process  and  is  essentially  an  extention  of  the 
frequency-domain  analysis  called  the  Method  of  Lines. 

The  Method  of  Lines  technique  was  developed  by  mathematicians  in  order  to 
solve  partial  differential  equations[27].  The  Method  of  Lines  is  simple  in  concept:  for 
a  given  system  of  partial  differential  equations,  all  but  one  of  the  independent 
variables  are  discretized  to  obtain  a  system  of  o.  dinary  differential  equations  so  that 
the  whole  space  is  represented  by  a  number  of  lines.  This  semi-analytical  procedure 
is  apparently  very  useful  in  the  calculation  of  planar  transmission  structures.  This  is 
because  these  structures  consist  of  regions,  which  are  homogeneous  in  one  direction 
as  shown  in  Fig.  1.3.  Moreover,  this  method  has  no  problem  with  so  called  "relative 
convergence"  phenomenon  sometimes  encountered  in  mode-matching  and  Galerkin's 


method. 
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The  application  of  Method  of  Lines  in  microwave  field  analysis  was  first 
introduced  in  [28]  for  the  calculation  of  the  propagation  constant  and  the  characteristic 
impedance  of  microstrip  lines.  The  method  has  also  been  extended  to  three 
dimensional  problems  such  as  the  calculation  of  the  resonant  frequencies  of  an 
arbitrarily  shaped  planar  resonator  [29].  There  have  also  been  several  publications  in 
which  the  Method  of  Lines  was  modified  for  the  analysis  of  a  quasi-planar  structure 
or  for  better  performance  [30].  Those  papers,  however,  have  always  used  the 
Method  of  Lines  to  solve  time-harmonic  Maxwell's  equations. 

This  study  describes  a  new  time-domain  method  termed  Time-Domain  Method 
of  Lines  (TDML).  The  goals  of  this  dissertation  are  as  follows. 

(1)  Development  of  the  procedure  and  formulation  of  the  new  method. 

(2)  Verification  of  the  method  by  solving  several  simple  two-dimensional  problems 
and  comparing  the  results  with  those  obtained  by  other  methods. 

(3)  Application  of  the  new  method  to  three-dimensional  structures  to  obtain  useful 
design  data. 

1.2  Dissertation  Organization 

Following  Chapter  1,  Chapter  2  presents  the  basic  discrete  mathematics 
related  to  the  Method  of  Lines  -  boundary  condition,  discretization  scheme,  difference 
operators  and  diagonalization  matrix.  In  Chapter  3,  the  analysis  of  the  two- 
dimensional  structure  with  a  uniform  interface  boundary  condition  is  formulated  and 
the  behavior  of  an  initial  pulse  inside  the  partially  filled  rectangular  waveguide  is 
investigated  using  the  present  method.  The  cutoff  frequency  spectrum  of  the  structure 
is  obtained  via  the  Fourier  transform  of  the  time-domain  data.  In  Chanter  4,  the 
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procedure  developed  in  Chapter  3  is  modified  to  deal  with  problems  where  the 
metalization  exists  on  the  aubsirate  surface.  In  this  situation,  a  non-uniform  interface 
boundary  condition  results.  The  procedure  is  used  to  calculate  the  cutoff  frequency 
spectrum  of  a  finned  waveguide  by  the  Fourier  transform  of  the  time-domain  data.  In 
Chapter  5,  the  formulation  of  the  Time-Domain  Method  of  Lines  for  the 
characterization  of  the  wave  propagation  and  scattering  property  in  the  three- 
dimensional  problems  is  described  in  detail.  The  results  of  the  time-domain  data  and 
the  derived  frequency-domain  characteristics  for  the  uniform  microstrip  line  and  its 
discontinuities  (step-in-width,  open-end  and  gap  discontinuities)  are  presented. 
Chapter  6  summarizes  the  contributions  of  this  study  and  proposes  some  related 
problems  for  future  research. 


Chapter  2 
Preliminary 


2.1.  Outline  of  the  Proposed  Method 

As  pointed  out  in  Chapter  1 ,  the  proposed  method  originates  from  the  fact  that 
most  of  the  discontinuites  appearing  in  the  planar  transmission  line  structures  are 
located  on  the  substrate  surface,  and  the  space  below  and  above  this  surface  is 
uniform  and  homogeneous.  One  wants  to  use  this  structural  feature  and  wishes  to 
solve  the  problem  by  discretizing  the  space  only  in  the  two-dimensional  surface  on  the 
substrate  where  the  discontinuity  is  located.  This  is  possible  if  the  wave-scattering 
information  in  the  direction  perpendicular  to  the  substrate  surface  is  available 
analytically.  The  new  method  actually  incorporates  this  process  and  entails 
discretization  of  the  structure  by  a  number  of  lines  perpendicular  to  substrate  surface 
as  shown  in  Fig  2.1.  At  a  specified  time,  the  field  distribution  at  each  intersection  of 
these  lines  with  the  substrate  surface  is  calculated  by  Maxwell's  equations  discretized 
only  in  the  x  and  z  directions  parallel  to  the  substrate  surface.  The  field  information  in 
the  y  direction  is  obtained  analytically  at  each  point  and  time.  This  information  can  be 
found  from  thc  inverse  spatial  transform  of  the  solution  of  the  frequency-domain 
Helmholtz  equation  in  the  y  direction. 

One  may  wonder  as  to  what  is  happening  to  the  wave  scattering  phenomena 
that  are  taking  place  everywhere  in  the  waveguide,  not  only  on  the  substrate  surface. 
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Fig.  2. 1 .  The  discretization  scheme  of  the  method  of  lines.  The  three-dimensional 
space  is  discretized  by  a  collection  of  lines  parallel  to  the  y-axis. 

This  question  is  natural,  because  in  other  time-domain  methods  the  electromagnetic 
fields  at  one  mesh  point  interact  with  those  at  all  six  neighboring  mesh  points  in  the  x, 
y,  and  z  directions.  In  the  proposed  method,  the  fields  at  any  point  on  one 
discretization  line  do  not  appear  to  interact  with  those  on  a  similar  point  on  another 
line.  It  should  be  emphasized  that  this  is  not  the  case.  As  will  be  seen  in  the 
formulation,  a  spatial  transformation  is  introduced  by  which  the  field  as  a  function  of 
(discretized)  x  and  z  is  transformed  to  another  discretized  quantity  which  contains  the 
field  quantities  at  all  x  and  z  values.  The  analytical  information  in  the  y  direction  is 
then  applied  to  this  transformed  quantity.  Since  analytical  expressions  are  used  for 
the  field  variation  in  the  y  direction,  this  method  can  easily  handle  the  case  where  the 
top  wall  is  removed  and  the  structure  is  open  in  the  y  direction. 


2.2.  Discrete  Mathematics  Related  to  the  Method  of  Lines 

As  explained  above,  the  proposed  method  is  a  "semi-analytical  method"  for 
solving  the  time-dependent  Maxwell's  equation  in  a  planar  waveguide  structure  so 
that  the  discrete  expressions  as  well  as  the  analytical  expressions  need  to  be  used  in 
the  formulation.  Therefore,  it  is  important  to  become  farmilar  with  some  discrete 
mathematics  related  to  the  method  of  lines  before  the  formulation  of  the  method  is 
explained.  This  mathematics  will  be  used  to  derive  all  the  formulations  described  in 
later  chapters. 

2.2.1.  Lateral  Boundary  Conditions  and  Discretization  Scheme 

Let  <J>  be  a  field  component  of  interest.  The  boundary  condition  for  4>  is 
always  given  by  either  the  Dirichlet  or  the  Neumann  condition  because  the  circuit  to 
be  considered  in  this  study  is  always  enclosed  by  perfect  electric  walls.  Since  all  the 
space  variables  (x  and  z)  except  one  (y)  are  discretized  in  the  Method  of  Lines,  the 
original  continuous  space  is  approximated  by  a  number  of  lines  parallel  to  the  y-axis 
located  only  at  discretized  positions  in  x-z  plane,  as  shown  in  Fig.  2.1.  The  locations 
of  these  discretized  points  are  determined  by  the  boundary  condition  of  the  field 
component,  0.  For  example,  if  a  boundary  is  known  to  be  a  Dirichlet-boundary  to  <J>, 
the  position  of  the  first  discretization  point  for  0  is  A£,  away  from  the  boundary  where 
A£  is  the  unit  discretization  length  in  ^  direction.  On  the  other  hand,  if  the  boundary 
is  a  Neumann-boundary,  the  first  discretization  position  is  A^/2  away  from  the 
boundary.  Figure  2.2  shows  the  difference  of  the  discretization  scheme  between  the 
Dirichlet-boundary  and  Neumann-boundary  .  One  thing  to  be  noticed  here  is  that  the 
two  set  of  the  discretization  points  are  shifted  by  A^/2.  This  discretization  scheme 


reduces  die  discretization  error  and  makes  the  second-order  difference  operator 
symmetric  [28]. 


4>  =  0  641/84  =  0 


Fig.  2.2.  Different  discretization  schemes  corresponding  to  the  boundary  conditions, 
(a)  Dirichlet-boundary,  (b)  Neumann-boundary. 

Figure  2.3  shows  the  outer  boundary  of  two  typical  structures  to  be 
considered  in  this  dissertation.  One  is  a  rectangular  waveguide  for  the  two- 
dimensional  problem  and  the  other  is  a  rectangular  waveguide  resonator  for  the  three- 
dimensional  problem.  The  waveguide  walls  consist  of  perfect  electric  walls  so  that 
the  lateral  boundary  condition  for  0  along  any  axis  is  always  given  by  either  the 
Dirichlet  or  Neumann  boundary  condition.  Since  there  always  exist  two  boundaries 
along  one  axis,  there  are  four  combinations  of  the  lateral  boundary  condition  in  one 
direction.  Figure  2.4  shows  the  diagrams  of  these  combinations  and  the  appropriate 
discretization  positions  for  each  case,  in  which  the  axis  is  assumed  to  be  x. 


Fig.  2.3.  The  outer  boundary  of  the  structures  considered  in  this  study. 


(a)  rectangular  waveguide,  (b)  rectangular  waveguide  resonator. 


Fig.  2.4.  Four  cases  of  lateral  boundary  condidon  combination. 


So  far,  only  the  one-dimensional  (along  the  x-axis)  discretization  scheme  has 
been  discussed.  However,  the  two-dimensional  discretization  (in  the  x-z  plane)  is 
required  to  solve  the  three-dimensional  problem.  In  that  case,  both  the  x-  and  the  z- 
direction  are  discretized  independently  according  to  the  discretization  scheme 
discribed  above.  A  typical  two-dimensional  discretization  example  is  shown  in  Fig. 
2.5  where  the  N-D  boundary  condition  is  assumed  along  both  the  x-direction  and  the 
z-direction. 


Fig.  2.5.  A  two-dimensional  discretization  example  where  N-D  boundary  condition 
is  assumed  along  the  x-  and  the  z-axis. 


2.2.2.  Matrix  Operators 


In  addition  to  the  discretization  of  the  space  according  to  the  boundary 
condition  of  the  Field,  the  derivatives  of  the  Field  with  respect  to  the  discretized 
variable  should  be  modiFied  to  a  discrete  form.  In  this  section,  the  difference 


operators  are  introduced  as  a  discrete  form  of  a  derivative,  and  the  relation  between 
them  will  be  derived. 
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2.2.2. 1.  First-Order  Difference  Operator 

After  the  discretization  of  the  space  using  the  discretization  scheme  discribed 
in  Section  2.2.1,  the  original  field  defined  in  the  continuous  whole  region  can  be 
approximated  by  a  collection  of  constituent  field  lines,  0j(y),  defined  only  at  the 

discretized  position  i.  In  the  cases  of  Fig.  2.4, 

<t>(x,y)  — >  [4)i(y)]1  =  02-  -  -  0n]1  (2.1) 


where  the  i-th  element  represents  the  field  of  the  i-th  line. 

Accordingly,  the  derivative  with  respect  to  the  discretized  variable,  x,  is 
replaced  by  the  difference  operation  of  two  consecutive  lines  using  the  central 
difference  scheme.  That  is, 

30(x,y) /3x --->  d0/9x  lj  =  [<J>j  -  0j.j]  /  Ax  fori=ltoN+l  (2.2) 

which  can  be  represented  in  terms  of  the  matrix. 

(2.3) 


A  d0(x,y) 

Ax  — =: - > 

dx 


1  1 
-1  1 


0 


0 


-1 1 


0o 

0i 


0N 

0N+1 


where  [Dx]  is  a  matrix  that  acts  as  a  difference  operator  with  respect  to  x. 


»  condition 
boundaryX 

Dirichlet 

Neumann 

left  boundary 

o 

o 

II 

o 

00  =  01 

right  boundary 

0N+1  =  0 

0N  =  0N+1 

Table  2.1.  Boundary  conditions  represented  by  the  field  values  near  the  boundary. 

The  specific  forms  of  the  difference  operator,  which  include  the  lateral 
boundary  conditions,  can  be  derived  by  imposing  the  appropriate  conditions  (see 
Table  2.1)  on  equation  (2.3).  The  resulting  first-order  difference  operators  are 
summarized  in  Table  2.2.  Although  the  operators  shown  in  Table  2.2  are  derived 
from  the  derivative  with  respect  to  x,  they  can  also  be  used  for  derivatives  with 
respect  to  the  other  variables.  From  table  2.2,  it  is  easy  to  find  the  relationship 
between  the  operators  which  will  be  used  in  the  later  formulation. 

[DN0]  = .  pDN]1  (2.4) 

[D^l  =  -  [Ddd]1  (2.5) 

where  the  superscripts  ND,DN,NN,and  DD  represent  the  specific  lateral  boundary 
conditions(left-right)  associated  with  the  difference  operator,  and  the  superscript  t 
represents  the  transpose  of  the  matrix. 

One  thing  to  be  noticed  here  is  that  the  values  of  the  derivative  obtained  from 
the  difference  operation  are  defined  only  at  the  mid-point  between  two  consecutive 


points  (i.e.  the  central  difference).  That  is,  the  actual  position  of  the  i-th  line  of  80  / 
dx  is  at  the  center  of  the  i-  and  (i+l)-th  positions  of  0  in  the  equation  (2.2). 


boundary 
conditions 
(left  -  right) 


elements  of  [d] 

di,i  = 

-  2  sin[(i-0.5)  n/ 
(2N+1)] 

(i  =  1 . .  N) 

di,i  = 

2  sin[(i-0.5)  rt  / 
(2N+1)] 

(i  =  1 . N) 

di+l,i 

=  2  sin[i  7t  / 

(2N+2)] 

(i  =  1. . ,  N) 

di-l,i 

=  2  sin[i  7t  / 

(2N+2)] 

(i  =  2, . .  N) 

Table  2.2.  First-order  difference  matrices  in  the  original  and  the  transformed  domain 
for  the  various  boundary  conditions. 


Therefore,  the  boundary  conditions  of  the  field  obtained  by  taking  a  difference 
operation  of  a  original  field,  (30  /  dx),  becomes  dual  condition  of  the  boundary 
condition  of  the  original  field,  <j>.  That  is,  if  the  boundary  condition  of  the  original 
field,  4>,  is  Neumann-Dirichlet  along  the  x  direction,  the  boundary  condition  of  its 
derivative  with  respect  to  x  becomes  Dirichlet-Neumann  along  the  x  direction. 

In  the  two-dimensional  discretization  of  the  three-dimensional  problem,  as 
shown  in  Fig.  2.5,  the  field  of  the  whole  region  is  considered  as  a  collection  of  Nx  * 
Nz  constituent  line  fields  which  is  represented  by  an  Nx  by  Nz  matrix  as  follows. 

$1.1  •  $l.Nz 

4>(x,y,z)  -->  [<t>i  /y)]  = 

$Nx.l  •  •  $Nx.N2 

where  the  <h,j  represents  the  field  of  the  (i,j)-th  line. 

For  the  first-order  derivative  of  with  respect  to  the  x-direction,  one  obtains 

3<t>(x,y,z)  90  $i+i,f$i,j 

dx  dx  A 

z  =  (j-l/2)Az  ,i=l . Nx  (2.7) 

or,  in  matrix  notation: 
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^l.Nz 

=  [DJ[0,jl 

^Nx.Mz 


(2.8) 


The  difference  matrix  [Dx]  depends  on  the  lateral  boundary  condition  for  <{>. 
As  shown  above,  it  has  the  same  form  of  operator  matrix  as  used  in  the  case  of  one¬ 
dimensional  discretization.  Here,  it  forms  the  difference  between  two  successive 
rows  of  matrix  [ty  j]. 

Analogously,  the  difference  operator  for  the  first-order  derivative  of  0  with 

respect  to  the  z-direction  should  form  the  difference  between  two  successive  columns 
of  the  matrix  [<5>i j]-  Thus,  the  difference  matrix  [DJ,  as  taken  from  the  Table  2.2  for 
the  N-D  boundary  condition  with  the  size  N  =  Nz,  will  operate  on  the  transpose  of  the 
matrix  fOjj). 


(2.9) 

(2.10) 


The  expression  in  equation  (2.10)  is  compatible  w  ith  the  difference  operators  for  the 
x-direction.  For  example. 


(2.11) 
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2. 2. 2. 2.  Second-Order  Difference  Operator 

Using  the  central  difference  scheme,  the  second-order  derivative  with  respect 
to  x  is  approximated  by 


0‘'4>(x,y)  >  32  4> 

dx2  d*2; 


0i+1-20i+<});.i 


Ax 


for  i  =  1  to  N 


(2.12) 


which  can  be  written  in  a  matrix  form. 


Ax 


2  cT4>(x,y) 
dx 


-> 


1  -2  1 
1  -2  1 


0 


0 


1  -2  1 


0o 

01 


0N 

0N+1 


=  -Pxx][0, 


(2.13) 


where  [Dxx]  represents  the  second-order  difference  operator  wth  respect  to  x.  The 
second-order  difference  operators  shown  in  Table  2.3  are  obtained  by  imposing  the 
boundary  conditions  of  the  field  ,  Table  2. 1 ,  to  equation  (2. 1 3). 

The  same  results  for  the  second-order  difference  operators  can  be  easily 
obtained  by  using  the  dual  boundary  condition  property  of  the  derivative  operation 
and  equations  (2.4)  and  (2.5);  that  is, 


I 


Ax  2  92<J>  /  9x2 


[d“]  [d^]  [oj  =  - 

_[d.1d!1N- 


[d“]1[d“]N  = 


fn00]  L 

-[Dx*J  <?>,] 


-MW 

-  M  [♦; 

-MK 


-(2.14a) 


(2.14b) 


(2.14c) 


(2. 14d) 


Although  the  results  shown  in  Table  2.3  are  derived  from  the  derivative  with  respect 
to  x,  they  can  be  used  for  the  derivatives  with  respective  to  the  other  variables. 

2.2.3.  Diagonalization  Transformation  Matrix 

Since  the  second-order  difference  operators  with  respect  to  x,  [Dxx],  are  real 

symmetric  matrices  as  shown  in  Table  2.3,  they  can  be  transformed  into  the  diagonal 
matrices  by  means  of  orthonormal  transformation  ma tries  [Tx]: 


[Tx00]1  [Dxxdd]  [Txdd]  =  diag  [dxxDDJ 
[Tx™]1  [Dxx™]  [TxNN]=diag  [dxxNN] 
[TxDN]t  [DxxDN]  f^DN]  =  diag  [dxxDN] 
[Tx^]1  [Dxx™]  [TXND]  =  diag  [dxxND] 


(2.15a) 
(2.15b) 
(2.15c) 
(2. 15d) 


rrxDD]t  [TxDD]  =  [TxNN]t  ITX^]  =  [TxDN]t  [J^N]  =  ^NDjt  ^NDj  =  INx 

(2.16) 


I 


boundary 

conditions 

left-right 

Dxx 

Tij 

N-D 

Of 

V  2  /  (N+0.5)  cos[(i-0.5)  (j-0.5)  n  / 
(N+0.5)] 

(i,j  =  l, . N) 

D-N 

V  2  /  (N+0.5)  sin[(i)  (j-0.5)  n  / 
(N+0.5)] 

(i,j  =  l, . N) 

D-D 

V  2/(N+l)sin[(i)(j)K  /(N+l)] 

(i,j  =  1 . N) 

N-N 

Bl 

V  2 / N cos[(i-0.5)  (j- DJt/N];j>l 

V  2  /  N  ;  j  =  1 

(i ,  j  =  1 . N) 

Table  2.3.  Second-order  difference  matrices  and  elements  of  the  transformation 
matrices. 


where  INx  represents  the  Nx  by  Nx  unit  matrix,  and  "diag"  means  a  diagonal  matrix 
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As  indicated  in  equations  (2.15a~d),  the  transformation  matrices  are  different 
depending  on  the  boundary  condition  of  the  second-order  difference  operator.  The 
formula  for  the  calculation  of  the  elements  of  the  transformation  matrices  are 
summarized  in  Table  2.3  (see  Apendix  1  for  the  derivation). 

Similarly,  the  first-order  difference  operators  can  be  (quasi-)diagonalized  by 
using  the  same  transformation  matrices  (see  Appendix  2  for  the  derivation) 


diag 


DD 

X 


(2.17a) 


r_DDi 

t 

C  nnI 

i_nni 

iT*  J 

[Dx 

T 

x 

— 

0 


(2.17b) 


=  diag 


DN 

X 


(2.17c) 


=  diag 


ND 

X 


(2.17d) 


Comparison  of  equations  (2.17a~d)  with  equations  (2.15a~d)  shows  that  the 
diagonalization  of  the  first-order  difference  matrix  requires  the  pre-multiplication  of 
the  transpose  of  the  transformation  matrix  for  the  dual  boundary  condition.  The 
forms  of  the  diagonalized  matrices  and  the  formula  for  the  calculation  of  the  diagonal 
elements  are  shown  in  Table  2.2. 

The  diagonal  elements  of  the  diagonalized  second-order  difference  matrices 
can  be  calculated  using  equations(2.14a~d)  and  (2.16).  That  is. 


diag  [dxxDD]  =  [TXDD]1  [DxxDDj  [^DDj 
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=  -  [TxDD]t  [DxLD]t  \TX™]\TX™]'  [Dxdd]  [Txdd] 

=  -  [dxD°]t  [dxDD]  (2.18a) 

diag  [dxxNN]  =  [TxNN3t  [D^NN]  [TxNN] 

=  -  [TxNN]t  [Tx^irTx00]1  PXNN]  [TXNN] 

=  -  [dxNN]t  [dxNNl  (2.18b) 

diag  [dxxDN]  =  [TxDN]t  [D^DN]  [T^N] 

=  -  [rxDN]£  PXDN]1  [TxND][TxND]t  PXDN]  [TXDN] 

=  -[dxDN]t(dxDN]  (2.18c) 

diag  [dxxK»]  =  [TxNDjt  [DxxND]  [TxND] 

=  -  [Tx^]1  [DxND]t  [TxDN][TxDN]t  [DXND]  [TxNDj 
=  -[dxNDlt[dxND]  (2.18d) 

2.3.  Summary 

In  this  chapter,  the  outline  of  the  proposed  method  was  described  briefly  and 
the  discrete  mathetics  related  to  the  method  of  lines  was  explained.  The  proposed 
method  makes  use  of  the  structural  feature  of  a  planar  circuit:  that  is,  the  spaces  above 
and  below  the  substrate  surface  are  uniform  and  homogeneous.  As  a  result,  the  space 
is  discretized  only  in  the  twodimensional  substrate  surface  and  an  analytical  solution 
is  used  for  the  normal  direction. 

The  descretization  scheme  used  in  the  method  of  lines  depends  on  the 
boundary  condition  of  the  field  of  interest.  This  shifted  descretization  scheme  reduces 
the  discretization  error  and  makes  the  formulation  symmetric.  Also,  the  discrete 
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forms  of  the  first-order  derivative  and  the  second-order  derivative  were  derived,  in 
which  the  side  wall  boundary  conditions  were  incorporated.  Finally,  the 
diagonalization  of  the  first-order  difference  matrix  and  the  second-order  difference 
matrix  were  explained  by  orthogonal  transformation  matrices. 


Chapter  3 

Time-Domain  Method  of  Lines  Applied  to  a  Partially  Filled 

Waveguide 


3.1.  Introduction 

Since  this  is  the  first  time  the  Time-Domain  Method  of  Lines  has  been  used  in 
the  solution  of  the  time-dependent  Maxwell’s  equations,  a  simple  problem  will  be 
used  as  a  test  case  to  check  the  validity  of  the  method.  The  simplest  possilble 
structure  to  analyze  is  the  homogeneously  filled  rectangular  waveguide  excited  by  an 
electric  field,  Ez,  infinite  in  length  and  uniform  in  the  z  (axial)  direction  so  that  a 
simple  two-dimensional  problem  results.  However,  the  formulation  to  be  described 
in  this  chapter  is  for  the  partially  filled  rectangular  waveguide  structure  shown  in  Fig. 
3.1.  This  is  appropriate  because  the  solution  to  the  homogeneously  filled  rectangular 


Fig.  3.1  Structure  of  a  partially  filled  rectangular  waveguide. 
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waveguide  can  be  obtained  as  a  special  case  of  the  partially  filled  waveguide  by  using 
the  same  dielectric  constant  for  region  I  and  region  II. 

3.2.  Formulation 

As  shown  in  Fig.  3.1,  the  problem  is  to  find  a  behavior  of  an  input  pulse  in  a 
partially  filled  rectangular  waveguide  excited  by  an  electric  field,  Ez,  infinite  in  length 

and  uniform  in  the  z  (axial)  direction.  The  structure  can  be  reduced  to  a  two- 
dimensional  one  because  of  the  input  condition.  The  formulation  for  such  a  problem 
is  simple,  yet  it  contains  all  the  essential  features  of  the  proposed  method.  This 
problem  corresponds  to  finding  the  cutoff  frequencies  of  various  TM  to  z  modes  in 
the  frequency  domain  [31].  Such  information  can  be  extracted  from  the  time-domain 
data  obtained  by  the  method  described  below. 

The  starting  point  of  formulation  is  the  time-dependent  Maxwell’s  equations 
expressed  in  scalar  notation;  that  is. 


dEz  /  dy  -  9Ey  /  3z  =  - 

p.  9HX  /  d  t 

(3.1a) 

3EX  /  dz  -  8EZ  /  3x  =  - 

p0Hy/3t 

(3.1b) 

dEy  /  3x  -  3EX  /3y  =  - 

■  p  3HZ  /  dt 

(3.1c) 

3Hz/3y-3Hy/az  = 

e(y)  3EX  /  dt 

(3. Id) 

8HX  /  dz  -  3HZ  /  3x  = 

e(y)  8Ey  /  3t 

(3.1e) 

dHy/dx-dHx/dy  = 

e(y)  dEz  /  dt 

(3.10 

Because  of  the  excitation,  only  Ez,  Hx  and  Hy  exist  and  d  /  dz  =  0.  Then, 


equations  (3.1a~f)  are  reduced  to 
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-  fi  aHx  /  at  =  BEZ  /  By 

(3.2a) 

-\i BUy/Bt  =  -BEz/Bx 

(3.2b) 

e(y)  3EZ  /  Bt  =  3Hy  /  3x  -  3HX  /  By 

'(3.2c) 

The  time-dependent  wave  equation  for  the  Ez  field  derived  from  Maxwell  s  equations 
(3.2a~c)  is 

32  Ez  /  3x2  +  32EZ  /  32y  -  M£( y)  32  Ez  /  32t  =  0  (3.3) 

Notice  that  the  side  wall  boundary  condition  for  Ez  and  Hx  is  Dirichlet- 
Dirichlet  and  that  of  Hy  is  Neumann-Neumann.  Figure  3.2  shows  the  cross-sectional 


Y 


Fig.  3.2.  Discretization  lines  for  Ez,  Hx,  and  Hy  incorporating  the  boundary 


conditions. 
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view  of  the  structure  shown  in  Fig.  3.1  with  the  appropriate  discretized  field  lines 
determined  by  the  rule  described  in  Section  2.2.1.  Now  the  time-domain  Maxwell's 
equations  and  the  wave  equation  are  discretized  in  the  x-direction  only. 


-^3[Hx]/3t  =  a[Ez]/ay  (3.4a) 

-\l  3[Hyj  /  at  =  -  [Dx™]  [Fz]  /  Ax  (3.4b) 

e(y)  atEJ  /  at  =  pxNN]  [Hy]  /  Ax  -  a[Hx]  /  By  (3.4c) 

[DxxDD][Ez]  /  (Ax)2  +  a2[Ez]  /  ay2  -  |i£(y)  B^E^  /  Bi  2=  0  (3.5) 


where  [Dxdd],  [DXNN]  are  the  first-order  difference  operators  and  [DXXDD]  is  the 
second-order  difference  operator,  in  which  the  side  wall  boundary  condition  is 
incorporated  as  explained  in  Section  2.2.2.  The  variables  [Ez],  [Hx]  and  [Hy]  are 
column  vectors  with  elements  which  are  actually  functions  of  y  and  t.  The  i-th 
element  of  the  column  vector  represents  the  field  along  i-th  line  in  Fig.  3.2. 

If  the  solution  for  [Ez]  is  found  using  equation  (3.5),  the  other  field 
components  are  obtained  from  the  equations  (3.4a~c).  However,  it  is  veiy  difficult  to 
solve  the  wave  equation  (3.5)  directly  because  the  equation  represents  a  set  of  coupled 
partial  differential  equations  due  to  [Dxxdd].  Since  [DXXDD]  is  a  real  symmetric 
matrix,  there  exists  a  real  orthonormal  matrix  [TXDD]  that  transforms  [DXXDD]  into  a 
diagonal  matrix  [dxxDD]  as  explained  in  Section  2.2.3. 

Using  equations  (2.15a)  and  (2.17a  and  b),  the  coupled  partial  differential 


equations,  (3.4a~c)  and  (3.5),  can  now  be  transformed  into  a  diagonalized  form  as 
follows. 
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(3.6a) 


(3.6b) 


(3.6c) 


(3.7) 


where  -  represents  the  transformed  field  corresponding  to  the  real  field  without  ~. 
The  transformation  of  the  real  field  is  accompanied  by  the  diagonalization  of  the  first- 
order  and  the  second-order  difference  matrices. 


N*[Tf 

N 

‘[HJ 

t 

II 
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tDD 
.  x 

(3.8a) 

(3.8b) 

(3.8c) 


Notice  that  equation  (3.7)  represents  a  set  of  uncoupled  partial  differential  equations. 
That  is,  it  can  be  solved  independently  along  the  i-th  line.  A  typical  equation 
describing  the  transformed  field  along  the  i-th  line  is 


Ez.i  + 


,DD 

-^7  E~  -  H£(y) 


Ax 


(3.9) 


with  the  boundary  conditions. 


EZil<y=0)=  Ez>1<y=b)  =  0 
Ez.i<y=h')  =  Hz,i(y=h+) 


3  1 

(3.10a) 

(3.10b) 


Using  the  separation  of  variables  technique,  one  can  obtain  the  solution  for  the  i-th 
line 


Ezi(y,0  = 

f  £n  (Ani  cos  conit  +  Bni  sin  tOnit)  sin  Klmfb-v) ,  for  region  I 


L  In  (Ani  cos  C0nit  +  Bni  sin  conit)  (sin 


for  region  II  (3.11) 


where  Klni,  K2ni  and  coni  are  determined  by  the  characteristic  transcendental 
equations  derived  from  the  interface  boundary  condition  (3.10b). 


Klni  cos  Klnjd  sin  K2nih  +  K2ni  sin  Klnid  cos  K2nih  —  0  (3.12a) 

o>ni2  =  [(Klni)2  -  dxxiDD  /  (Ax)2]  /  m 

=  [(K2ni)  2  -  dxxlDD  /  (Ax)2]  /  (3.12b) 


Notice  that  for  a  given  n  in  equation  (3.1 1),  the  underlined  part  constitutes  the  n-th 
mode  on  the  i-th  line,  ^zni  ^  . 


EZm(y)  = 
r  sin  Klni(b-y), 

L  (sin  Klnid  /  sin  K2nih)  sin  K2niy  , 


in  region  I 

in  region  II  (3.13) 


From  this  point  on,  there  are  basically  two  approaches.  First,  by  knowing  the 
initial  conditions  for  Ez,  one  can  find  Anj  using  the  orthogonal  property  of  the  modes 


3  2 


[32]. 


Ani  = 


I. 


Ezl(t=0)  ^oE(y)  Ezni<y)  dy 


I 


Ezn.<y)  ^oE(y)  dy 


(3.14) 


One  way  to  determine  Bnj  is  to  use  the  causality  condition  for  Hx  and  Hy, 
which  states 


H^i[y>t-‘  2)  2  ^  °  (3.15) 

where  At  is  the  unit  time  interval  in  the  sampling  process.  The  expressions  for  Hxi 
and  Hyi  can  be  found  by  integrating  equations  (3.6a  and  b)  with  respect  to  t. 


sin  co  -t 

n  i  .a  ni 


C0„ 


cos  Kni  y 


(3.16a) 


Hy  i  (t)  =  * 


p0Ax 


(3.16b) 
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,DD 

dx.» 

(IqAx 


K„i 


-  B. 


COS  CO  -t 


—  +  A. 


sin  conit 


to 

m 


CO  , 

n>  / 


sin  Kn)y 


Imposing  the  causality  condition,  equation  (3.15),  to  equations  (3.16a  and  b),  one 
can  obtain 


Bni  —  -Ani  tan  (coni  At  /  2) 


(3.17) 


After  Anj[  and  Bni  are  determined,  the  solution  at  any  point  and  time  can  be 
found  by  using  equation  (3.11).  Then,  the  real  field  can  be  obtained  by  the  inverse 
transform. 


E,(y,t)' 


tdd 

A  X 


Ez(y,t)] 


(3.18) 


Second,  an  alternative  method  is  the  application  of  a  time-stepping  procedure, 
where  equations  (3.6a~c)  are  discretized  in  time: 

H~(y,N  +  j)  =  H^,<y,N  -  h  -  —  A  E~(y,N) 

^0  ay  (3.19a) 

—  At  d00  _ 

H“(y,N  +  i)  =  H~(y,N  -  h  +  —  -2-  Ez<1<y,N) 

Ax  (3.19b) 

A  d00  At 

E7;(y,N+ 1 )  =  E~(y,N)  -  — - —  H~(y,N  +  j)  -  — - H^,<y,N  +  I) 

e(y)  Ax  z  £( y)dy  z 


(3.19c) 
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From  the  initial  condition  of  Ez,  one  can  find  Ani  at  time  t  =  0  using  equation 
(3.14).  We  denote  this  quantity  as  AniN  with  N  =  0.  The  solution  (3.11)  at  the  sime 
step  N  can  now  be  expressed  in  the  form 

Ezi  (y)  = 

f  Xn  AniN  sin  Klni(b-y), 

L  Xn  AniN  (sin  Klnia  /  sin  K2nih)  sin  K2niy, 

Similarly,  for  Hxi  and  Hyj  at  the  time  step  (N  +  1/2), 

JnTT/2. 

hx1  (y)  = 


r  In  BmN+1/2  cos  Klni(b-y), 

in  region  I 

L  Xn  BniN+1'-  (cos  Klnid  /  cos  K2njh)  cos  K2niy, 

in  region  II 

(3.20b) 

X 

^  ’ 

1 

II 

rXn  CmN+1/2  sinKlni(b-y), 

in  region  I 

L  Xn  CniN+1/2  (sin  Klnid  /  sin  K2nih)  sin  K2niy, 

in  region  II 

(3.20c) 

Substituting  the  equations  for  the  Fields,  (3.20a~c),  into  equations  (3.19a~c), 
a  leap-frog  type  iteration  scheme  [IF1,  can  be  implemented  to  calculate  the  coefficients 
of  modal  fields  on  the  i-th  line  at  any  time  step  N. 


in  region  I 

in  region  II  (3.20a) 


N+l/2  N-l/2 

n  i 


At  .N 
+  Kni  Anj 


(3.21a) 


N+l/2 

1 

At  DD  .  N 

ni 

=  cni  + 

—  d*,i  Anj 

^0 

Ax 

(3.21b) 

N+l  _ 

.  N  1 

At 

, NN  _N+l/2 

At  _  _  N+l/2 

- 

ni  — 

■^ni 

dx.iCni  - 

- K„.iBni 

e(y) 

Ax 

£(y) 

(3.21c) 

where  Knj  becomes  either  K 1  ni  or  K2nj  depending  on  the  region.  After  the 

coefficients  at  the  N-th  time  step  are  determined,  the  transformed  field  can  be  obtained 

by  the  equations  (3.20a~3).  The  real  field  at  the  N-th  time  step  can  be  obtained  by 

rg-]N 

invoking  the  inverse  transformation  as  described  above  to  I  A  . 

After  the  history  of  pulse  scattering  is  obtained  in  the  time-domain,  the  cutoff 
frequency  spectrum  of  the  structure  can  be  extracted  via  the  Fourier  transform  of  the 
time-domain  data.  If  the  time-domain  data  consists  of  N  data  points  with  unit 
sampling  time  At,  the  resolution  of  the  frequency  domain  information  is  Af  =  1  /  NAt 
,  and  the  maximum  frequency  range  that  can  be  obtained  is  NAf/2. 

For  comparison,  the  cutoff  spectrum  of  the  same  structure  can  also  be 
calculated  analytically  by  using  the  transverse  resonance  technique.  For  TM  to  z 
mode  (LSE  mode),  the  eigenvalue  equation  is  given  by  [32] 


p  tan  (qh)  =  -q  tan  (pd)  (3.22a) 

p2  =  £j  o^iiqEq  -  (n7t/a)2  (3.22b) 

q2  =  £2  co2(J.q£q  -  (nrt/a)2,  forn  =  l, .  (3.22c) 


3.3.  Results  and  Discussion 

3.3.1.  Homogeneously  Filled  Rectangular  Waveguide 
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The  structure  considered  has  dimensions:  a  =  2  [cm],  b  =  1  [cm],  ej  =  3; 

£2  =  3.  Figure  3.3  depicts  the  Ez  field  distributions  in  the  waveguide  cross-section 
(x-y  plane)  at  each  time  step  after  a  pulsed  Ez  excitation  is  imposed  at  t  =  0  at  the 
center  of  the  cross  section.  Figure  3.4  shows  the  spectrum  of  the  time  signal  of  Ez 
where  the  waveguide  cutoff  frequencies  are  represented  by  the  peaks.  The  cutoff 
frequencies  obtained  by  present  method  differ  by  less  than  one  percent  from  the 
analytical  values. 

3.3.2.  Partially  Filled  Rectangular  Waveguide 

The  structure  is  ;  a  =  2  [cm],  b  =  1  [cm],  h  =  0.3  [cm],  £j  -  ^  E2  =  3. 

Figure  3.5  depicts  the  Ez  field  distributions  in  the  waveguide  cross  section  (x-y  plane) 
after  a  pulsed  Ez  excitation  is  imposed  at  t  =  0  at  the  center  of  the  cross  section. 
Figure  3.6  shows  the  spectrum  of  the  time  signal  of  Ez  where  the  waveguide  cutoff 
frequencies  are  represented  by  the  peaks.  The  cutoff  frequencies  obtained  differ  by 
less  than  one  percent  from  the  analytical  values. 

It  should  be  noted  that  the  two  approaches  ( the  modal  expansion  and  the  time¬ 
stepping  )  described  above  should  be  equivalent  in  principle.  However,  each  has 
advantages  and  disadvantages.  For  instance,  if  one  deals  with  a  time-dependent 
excitation,  the  time-stepping  method  would  be  simpler  to  implement.  Otherwise,  the 
modal  expansion  method  is  more  efficient  as  long  as  only  the  results  at  a  particular 
time  are  of  interest.  However,  if  any  frequency-domain  information  is  needed,  the 
time  history  needs  to  be  found.  The  first  method  needs  to  be  used  at  many  time 
instances.  The  time-stepping  method  automatically  generates  the  time  history. 
Hence,  the  amount  of  computation  would  be  about  the  same.  In  many  cases,  there 


Fig.  3.3  The  Ez-pulse  scattering  in  a  homogeneously  filled  rectangular  waveguide 
(a=2,  b=l  cm,  e  =3).  (a)  t=0,  (b)  t=20,  (c)  t=30  and  (d)  t=80  p-sec. 


Fig.  3.4  The  cutoff  frequency  spectrum  of  the  homogeneously  filled  waveguide 


when  the  initial  input  is  placed  at  the  center  of  the  waveguide. 
(a)TMll,  (b)  TM31  and  (c)TM51  modes 


Fig.  3.5  The  Ez-pulse  scattering  in  a  partially  filled  rectangular  waveguide 


are  ambiguities  in  finding  the  initial  condition  for  the  time  derivative.  In  such 
instances,  it  may  be  simpler  to  use  time  stepping  to  generate  the  time  history 
required.  Finally,  it  is  also  possible  to  switch  from  one  method  to  the  other . 

3.4.  Conclusions 

A  new  time-domain  technique  is  presented  in  which  an  analytical  process  is 
incoporated  along  one  of  the  spatial  dimensions  so  that  the  dimensions  of  the  problem 
are  effectively  reduced  by  one.  The  present  approach  can  be  used  to  calculate  the 
cutoff  frequency  of  other  planar  transmission  structures  and  can  be  extended  to 
propagation  problems.  It  has  a  number  of  potential  advantages  over  many  other  time 
domain  methods.  First,  the  method  is  believed  to  be  efficient  since  much  analytical 
processing  is  used.  Second,  the  problem  can  handle  open  boundaries  in  the  vertical 
(y)  direction  because  of  the  analytical  solutions  used  in  y. 


Chapter  4 

Time-Domain  Method  of  Lines  Applied  to  Finned 

Waveguide 


4.1.  Introduction 

A  new  time-domain  method  was  studied  for  a  partially  filled  rectangular 
waveguide  in  the  Chapter  3  in  which  the  interface  boundary  condition  is  not  changed 
along  one  of  the  spatial  dimensions  (transformed  direction).  Therefore,  the  problem 
could  be  reduced  into  a  one-dimensional  boundary  value  problem  which  requires  only 
one-dimensional  discretization. 

In  the  previous  chapter,  the  feasibility  of  a  new  time-domain  method  was 
shown  by  calculating  the  cutoff  frequencies  of  the  partially  filled  waveguide  problem 
accurately.  In  this  chapter,  the  method  will  be  extended  to  a  structure  which  contains 
a  metallic  strip  at  the  interface  boundary.  This  extension  also  makes  use  of  the 
Method  of  Lines  and  an  analytical  process  can  be  incorporated  along  one  of  the  spatial 
dimensions.  However,  because  of  the  non-uniform  boundary  condition  in  the 
transformed  direction,  each  line  is  not  independent  but  related  to  each  other  in  order  to 
match  the  required  interface  boundary  conditions.  Consequently,  instead  of 
considering  one-dimensional  eigenmodes  for  each  line  as  in  the  case  of  the  uniform 
interface  boundary  condition  problem,  two-dimensional  eigenmodes  should  be 
considered  in  this  structure  [33]. 


40 


4.2.  Formulation 

As  shown  in  Fig.  4.1,  the  problem  is  to  find  the  time-domain  behavior  of  an 
input  pulse  in  a  uniform  rectangular  waveguide  excited  by  an  electric  field  infinite  in 
length  and  uniform  in  the  z-direction.  The  structure  can  be  reduced  to  a  two- 
dimensional  one  because  of  the  input  condition.  This  problem  corresponds  to  finding 
the  time-domain  behavior  of  the  pulsed  input  and  the  cutoff  spectrum  of  the  given 
structure.  The  formulation  is  very  similar  to  the  partially  filled  waveguide  problem 
except  for  the  treatment  of  the  interface  boundary  condition. 


Fig.  4.1  Structure  of  a  finned  rectangular  waveguide. 


Because  of  the  excitation,  only  Ez,  Hx  and  Hy  exist  and  d  /  dz  =  0.  Also, 
notice  that  the  side  wall  boundary  condition  of  the  fields  Ez,  Hx  is  Dirichlet-Dirichlet 
and  that  of  Hy  is  Neumann-Neumann.  Since  the  structure  is  uniform  in  the  space 
above  and  below  the  y  =  h  axis,  the  time-domain  equations  are  discretized  in  the  x- 
direction  only.  Figure  4.2  shows  the  cross-sectional  view  of  the  finned  rectangular 


waveguide  with  the  appropriate  discretized  field  lines  incoporating  the  boundary 
conditions.  Now,  Maxwell's  equations  in  (3.1  a— f)  are  reduced  and  discretized  to 


h  2 


Y 


Fig.  4.2.  The  cross-sectional  structure  of  a  finned  rectangular  waveguide. 


■^id[Hx]/dt=d[Ez]/dy  (4.1a) 

-p.  d[Hy]  /  8t  =  [DxDD]  [EJ  /  Ax  (4. 1  b) 

e(y)  a[Ez]  /  3t  =  [DXNN]  [Hy ]  /  Ax  -  3[Hx]  /  9y  (4.1c) 

The  governing  equation  for  Ez  is 

PxxDD][Ezl  /  (Ax)2  +  a^EJ  /  ay2  -  q£(y)  d2[Ez]  /  9t2  =  0  (4.2) 


with  appropriate  boundary  conditions  ; 


[Ez(t,y=h+)]  =  [Ez(t,y=h-)],  for  all  t  and  i 


(4.3a) 


[Hx(t,y=h-)]  -  [Hx(t,y=h+)]  =  [Jz(t,y=h)],  for  all  t  and  i 
[Ez(t,y=h)]  =  0,  for  all  t  and  i  on  M 
[Ez(t,y=0)]  =  [Ez(t,y=b)]  =  0,  for  all  t  and  i 


(4.3b) 

(4.3c) 

(4.3d) 


4  3 


where  [DXDD],  [D;*^],  [DXXDD]  are  difference  operators  in  which  the  side  wall 
boundary  condition  is  incorporated  (see  Tables  2.2  and  2.3).  The  variables  [Ez], 
[Hx]  and  [Hy]  are  column  vectors  representing  the  fields  along  each  line  and  are 
functions  of  y  and  t.  Notice  that  equation  (4.2)  is  a  system  of  a  coupled  partial 
differential  equations. 

Since  [D^00]  is  a  real  symmetric  matrix,  there  exists  a  real  orthogonal  matrix 
[TxDD]  that  transforms  [DXXDD]  into  a  diagonal  matrix  [dxxDD],  We  can  now 
transform  the  real  field  [EJ  into  a  transformed  field  [EJ  =  Hx00]1  [EJ  where  the 
superscript  t  stands  for  transpose.  The  transform  of  equation  (4.2)  is 


Md“]N  +  *  [Ej-neW®  [Ez]  =  [0] 

Ax  dy  dt 

(4.4) 

with  the  transformed  boundary  conditions  of  (4.3a~d) 

E^t,  y-h  )  -  E^t,  y-h  )  ,  fQr  j 

(4.5a) 

[«*«■  y-h*)]  -  find,  y=h")]  =  -[j'A  y-h)].  foraU ,  Md  i 

(4.5b) 

[^t.  y-h)]  =  unknown  .  for  a]{  <  and  i  on  M 

(4.5c) 

[E^y=0)]=[E*y=b)]  =  [0]>foralJIandi 

(4.5d) 
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Notice  that  equation  (4.4)  is  a  set  of  uncoupled  partial  differential  equations 
related  by  non-uniform  boundary  conditions  (4.3a~c).  For  the  i-th  line,  equations 


(4.4)  -  (4.5c)  become 

iCD  2  -x2 

-^Ez,i  +  ^TEZti-ne(y)^TEz,i=0 

Ax’  ^  3‘ 

(4.6) 

E^i(t,  y=h+)  =  E~(t,  y=h ) 

(4.7a) 

H^it,  y=h+)  -  y=h)  =  -  y=h) 

(4.7b) 

Ez  i(t,  y=h)  =  unknown 

(4.7c) 

yt,  y=0)  =  E~(t,  y=b)  =  0 

(4.7d) 

By  means  of  the  separation  of  variables  technique,  the  intermediate  solution  on  the  i- 
th  line  for  the  n-th  mode  which  satisfies  the  transformed  boundary  condition,  (4.5a 
and  d),  is 

E~(t,y)  = 

(Ani  cos  conit  +  Bni  sin  (Unit)  sin  Klni(b-y),  in  region  I 

(Ani  cos  conjt  +  Bni  sin  coniO  (sin  Klnid  /  sin  K2nih)  sin  K2niy,  in  region  II 

(4.8a) 

where  Klni,  K2ni,  and  coni  are  related  by 

Klni2  =  coni2  p£]  +  dxxiDD  /  (Ax)2 
K2ni2  =  o>ni2  ^i£2  +  dxxiDD  /  (Ax)2 


(4.8b) 


The  modal  current  [Jzn]  obtained  from  the  discontinuity  of  the  corresponding  modal 
magnetic  field  [Hxn]  at  y  =  h  by  using  (4.1a)  is 


J„.„,<t,y=h)  =  ■  ^  K2ni  cos  K2nih  +  K1  nicos  K1  nfJ 

=  Y^i(t.y=h)  £Wt>y=h) 


(4.9a) 


where  fni(t)  =  (Ani  sin  (Unit  /  CDni  -  Bni  cos  C0nit  /  CDni) 


(4.9b) 


1  1  Anisin(onit-  Bnicosci)nit 

p  coni  Anicosconit-Bnisino)nit 


(KlnicotKlnid-fK2nicotK2nih) 


(4.9c) 

Now,  the  final  boundary  condition,  (4.3b),  is  applied  in  the  real  field  domain. 


(4.10) 


where  r¥n(hy=h)]  is  the  diagonal  admitance  matrix  with  diagonal  elements  given  by 
equation  (4.9c).  Because  Jzni  =  0  if  i  is  on  the  non-metallic  portion  (D)  and  Ezni  =  0 
if  i  is  on  the  metallic  portion  (M),  the  characteristic  matrix  equation  can  be 
derived. from  equation  (4.10). 


(4.11) 


where  "D”  represents  the  non-metalization  portion  of  the  substrate  surface  as  shown 
in  Fig.  4.2  and  the  subscript  "red"  represents  the  reduced  matrix  obtained  from  the 
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full  matrix  in  {  }  by  removing  the  rows  and  columns  corresponding  to  the 
metalization  portion. 

In  order  for  the  n-th  eigenmode,  [Ezn],  to  satisfy  (4.11)  at  all  times  with  a 
non-trivial  solution,  the  time -dependent  solution  for  each  line  should  have  the  same 
functional  form  except  for  a  constant  factor,  and  the  determinant  of  the  reduced  matrix 
should  be  zero.  The  final  solution  of  the  [Ezn]  is 

E^.-O-y)  = 

[  (An  cos  C0nt  +  Bn  sin  (Ont)  Cni  sin  Klni(b-y),  in  region  I 
L  (An  cos  tOnt  +  Bn  sin  (Ont)  Cni  (sin  Klnjd  /  sin  K2nih)  sin  K2niy,  in  region  II 

(4.12) 

where  (On  is  n-th  eigenvalue  of  the  characteristic  equation  given  by 


(4.13) 


and  the  coefficients,  Cni,  can  be  derived  from  the  eigenvector,  [Ezn],  of  (4.11) 
corresponding  to  the  n-th  eigenvalue.  The  procedure  described  above  is  very  similar 
to  the  frequency-domain  method  of  lines  [28]. 

Finally,  any  real  field  can  be  expanded  in  terms  of  the  eigenmodes  in  the 
transformed  domain 

ETi(t-y)  = 

f  Xn  (An  cos  tont  +  Bn  sin  o>nt)  Cni  sin  Klm(b-v).  in  region  I 

L  In  (An  cos  <ont  +  Bn  sin  tont)  Cni  (sinKImd  /  sin  K2njh)  sin  K2njy , 


in  region  II 


(4.14) 
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Notice  that  for  a  given  n  in  equation  (4.14),  the  underlined  part  constitutes  the  n-th 
mode  of  the  i-th  line 


Ei.ni(y)  = 

|-SinKl„,(b-y)f  inregionI 


sin  Kln.d 

i  — 7-sinK2my 

l_sinK2nih  niJ 


,  in  region  II 


(4.15a) 

(4.15b) 


In  equation  (4.14),  An  can  be  found  from  the  initial  condition  of 
Ez(t-0,y)]  -  [tx  [E/t-O.y)]  orthogonality  property 


A 


n 


Ez(t=0,y) 


Ezn(y)]dy 


[Ezn(y)  ]  [Ezn(y)]  dy 


(4.16) 


Also,  the  causality  condition  assumed  by  time  iteration  method  [34]  can  be  used  to 
find  Bn 

Bn  =  An  tan(o)n  At  /  2),  where  At  <  Ax  /  V  2  Cjnax  (4.17) 

where  Cmax  is  the  maximum  wave  phase  velocity  within  the  structure.  The  real  field 
at  any  time  can  be  obtained  by  invoking  the  inverse  transformation  to  the  transformed 
Held, 


(4.18) 
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4.3.  Results  and  Discussion 

The  structure  of  the  finned  waveguide  considered  is  a  =  1  [cm],  b  =  2  [cm],  h 
=  0.5  [cm],  M  =  0.25  [cm],  ei  =  3,  and  £2  =  3.  To  check  the  validity  of  the  method, 
the  eigenfrequency  of  the  finned  rectangular  waveguide  for  the  dominant  TE  mode 
was  found  and  compared  with  Hoefer's  result  [24]  in  Fig.  4.3.  Good  agreement  was 


with  b/a=2,  h=l. 

( light  squares  ;  present  method,  dark  squares  ;  Hoefer  [24] ) 

obtained.  For  further  confirmation,  the  Ez  field  distributions  in  the  homogeneously 
filled  rectangular  waveguide  are  calculated  by  the  present  two-dimensional  method, 
and  the  result  is  compared  with  that  obtained  by  the  one-dimensional  method 
described  in  the  Chapter  3.  Figure  4.4  shows  the  lesult.  Even  though  the  ripple 
seems  to  be  higher  than  that  of  the  one-dimensionai  analysis,  the  result  shows  that  the 
wave  scattering  characteristics  can  be  analyzed  by  the  two-dimensional  time-domain 
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method.  The  pulse  propagation  and  scattering  in  the  finned  rectangular  waveguide  is 
dipicted  in  Fig.  4.5. 


m  i“i 


Fig.  4.4.  Pictures  of  a  pulse  scattering  in  a  partially  filled  rectangular  waveguide 
with  a=2,  b=l,  h=.3,  £1=1,  £2=3.  ([I] ;  one-dimensional  method, 

[H] ;  two-dimensional  method) 

(a)  t=0,  (b)  t=20,  (c)  t=40,  and  (d)  t=60  [p-sec] 

4.4.  Conclusions 

In  this  chapter,  it  was  shown  that  the  two-dimensional  Time-Domain  Method 
of  Lines  can  be  used  to  analyze  planar  transmission  structures  which  have  a  non- 
uniform  boundary  condition  at  the  interface  boundary  along  the  transformed  direction. 
This  is  accomplished  by  the  two-dimensional  eigenmode  expansion  concept  instead  of 
the  one-dimensional  method  used  in  the  case  of  the  uniform  boundary  problem. 


Chapter  5 

Characterization  of  Uniform  Microstrip  Line  and  its 
Discontinuities  using  Time-Domain  Method  of  Lines 


5.1.  Introduction 

In  the  previous  chapters,  the  TDML  has  been  applied  to  the  two-dimensional 
problems  to  obtain  the  scattering  data  of  the  input  pulse  in  the  time  domain.  The 
cutoff  frequency  spectrums  of  planar  transmission  lines  were  extracted  from  the  time- 
domain  data  via  the  Fourier  transform.  Although  the  cutoff  frequency  is  an  important 
parameter  characterizing  a  waveguide  structure,  the  propagation  characteristics  are 
even  more  important  to  the  microwave  circuit  designer.  In  this  chapter,  it  will  be 
shown  that  the  TDML  can  be  extended  to  the  three-dimensional  propagation 
problems.  The  time-domain  behavior  of  the  input  pulse  in  a  planar  transmission  line 
circuit  can  be  obtained  by  TDML.  From  the  time-domain  data,  the  frequency  domain 
characteristics  for  a  wide  range  of  frequencies  can  be  found  via  the  Fourier  transform. 
The  procedure  and  the  formulation  for  the  characterization  of  the  three-dimensional 
structure  is  presented  [35].  The  characteristics  of  a  uniform  microstrip  line  and  its 
discontinuities  (step-in-width,  open-end  and  gap  discontinuities)  are  presented  and 
compared  with  available  published  data. 
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_ 1 _ i 

input  field 


Fig.  5. 1  A  general  planar  transmission  line  with  discontinuity  in 
a  shielded  rectangular  waveguide 

5.2.  Formulation 

Let  us  consider  a  general  planar  transmission  line  with  a  discontinuity,  as 
shown  in  Fig.  5.1.  In  order  to  apply  the  TDML,  electric  or  magnetic  walls  need  to  be 
placed  at  both  ends  of  the  waveguide  so  that  a  resonant  structure  will  result.  Figure 
5.2  shows  the  top  view  of  the  structure  with  discretization  points  for  the  y  component 
of  the  electric  Field,  Ey.  It  is  assumed  that  the  structure  has  spatial  symmetry  in  x- 
direction  so  that  the  problem  domain  can  be  reduced  by  a  factor  of  two. 

5.2.1.  Discretization 

As  shown  in  Fig.  5.2,  the  structure  is  discretized  by  Field  lines  which  are 
properly  placed  to  satisfy  the  boundary  conditions  on  the  side  walls  and  the  end 


53 


side  wall  A  ;  M-wall 


Z  c  Nz  ....  2  j=l 


i=l 

2 


Nx 


a 


X 


Fig.  5.2  Top  view  of  the  structure  shown  in  Fig.5.1  with  the  proper 


discretization  points  for  Ey  component 


walls.  For  example,  since  the  Neumann  boundary  condition  for  Ey  is  applied  on  the 
side  wall  A  and  the  Dirichlet  condition  on  the  side  wall  B,  the  Ey  points  are  located 
away  from  the  side  wall  A  by  one  half  of  Ax  and  away  from  side  wall  B  by  one  Ax. 
Also,  it  is  important  to  choose  the  Ey  nodes  as  the  outermost  nodes  on  the  strip  in 
order  to  reduce  the  error  in  the  calculation  of  the  field  distribution.  After  discretization 
of  the  space,  the  original  field  in  the  continuous  space  is  approximated  by  many  field 
lines  parallel  to  the  y-axis  at  the  discretization  points 

Ey(x,y,z,t)  —  >  [Hy^kCy.t)]1  =  [Ey,i,  Ey,2, ....  Ey.NxNzl1  (5.1) 

1  <i  <NX,  1  <j  <NZ. 


where  k  =  i  +  (j-1)  Nx, 


(5.2) 
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5.2.2.  Expansion  of  the  Input  Pulse 

In  the  structure  shown  in  Fig.  5.1  with  end  walls,  a  field  distribution  at  any 
time  t,  [Ey(y,t)],  can  be  expanded  by  modal  field  distributions,  [Eyn(y)]  as  follows; 

[Ey(y,t)]  =  £n  (An  cos  0>nt  +  Bn  sin  CDnt)  [Eyn(y)]  (5.3) 

where  [  ]  denotes  a  column  vector  whose  i-th  component  represents  the  field  value  at 
the  i-th  discretization  line  as  shown  in  equation  (5.1).  Also  the  subscript  n  denotes 
the  mode  number.  One  thing  to  be  mentioned  here  is  that  if  the  structure  and  the 
initial  input  support  any  static  charge  distribution,  the  DC-mode  (con  =  0)  should  be 
included  in  equation  (5.3)  to  obtain  the  correct  time-  and  frequency-domain  data. 

In  the  expansion  of  (5.3),  the  coefficients  An  can  be  determined  by  the  input 
field  distribution,  [Ey(y,  t=0)]  and  the  orthogonal  property  of  the  modal  fields 


An  = 


f  [Ey(y,t=0)jEy  n(y)]dy 

J  0 

f  [Ey,n(y)]tEy.n(y)]dy 

J  o 


(5.4) 


Also,  the  Bn  can  be  found  by  the  causality  condition  used  in  the  time  iteration 
method[  19,  34] 


Bn  =  An  tan  (o>n  At/ 2) 
where  At  <  min(Ax,Az)  /  Cmax  ^  3 


(5.5) 

(5.6) 


and  Cmax  is  the  maximum  wave  phase  velocity  within  the  structure. 


In  order  to  find  the  engenfrequencies  and  eigenmodes  of  the  structure,  a 
technique  used  in  frequency-domain  Method  of  Lines  [29]  is  borrowed.  In  this 
technique,  the  characteristic  equation  is  obtained  by  applying  the  metallic  boundary 
condition  to  the  metalization  on  the  dielectric  interface  boundary.  The  equation  is 
given  by  (see  Apendix  3  for  details) 


[N1 

\M 

red 

[PJ1 

[p-al 

[p] 

red 

[N1 

[p«U 

red 

[Gzz][Gzx] 

_[Gxz][Gxx]_ 

[H 


for  static  fields 


(5.7a) 


for  time-harmonic  fields  (5.7b) 


In  order  to  obtain  a  non-trivial  solution  of  equations  (5.7a  and  b), 

det  [G(co)]  red  =  0,  for  time-harmonic  fields  (5.8) 

where  the  subscript  "red"  represents  the  reduced  matrix  corresponding  to  the 
discretization  points  on  the  metalization. 

After  the  eigenfrequencies  are  found,  equations  (5.7a)  and  (5.7b)  can  be  used 
to  obtain  the  charge  distribution,  [p]1,  and  the  current  density  distributions, 
[[JzUJx]!1-  The  modal  field  distributions,  [Eyn(y)l,  can  be  derived  from  these 
quantities.  It  is  found  that  the  system  of  equations  (5.7a)  and  (5.7b)  need  to  be 
solved  by  the  QR  method  [40]  to  ensure  stability  of  the  results. 


5.2.3.  Time-Domain  Data  and  Frequency-Domain  Characteristics 


Once  the  procedures  described  in  Sections  5.2.1  and  5.2.2  are  carried  out,  all 
the  constants  in  equation  (5.3)  are  determined  so  that  the  field  value  at  any  point  and 
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time  can  be  calculated.  Therefore,  the  characteristics  of  a  circuit  can  be  obtained  by 
observing  the  propagation  and  scattering  of  the  input  pulse  in  the  given  circuit.  In  this 
chapter,  four  structures  will  be  tested,  namely,  the  uniform  microstrip  line  and  three 
commonly  encountered  discontinuities  (step-in-width,  open-end  and  gap). 


— 

A 

L 

0 

Fig.  5.3  Top  view  of  a  uniform  microstrip  line  with  two  observation  planes 
and  initial  input  pulse. 

5.2.3. 1.  Uniform  Microstrip  line 

Figure  5.3  shows  a  uniform  microstrip  line  with  two  observation  planes.  In 
this  case,  the  transfer  function  can  be  found  from 

expf-y(o))Ll  =  V-('*  ^=i  t  /  Vy(co,z=0)  (5.9) 

where  yfco)  =  a(co)+j(3(o))  =  propagation  constant  of  the  uniform  microstrip  line. 
Vy«o,z=0)  =  Fourier  transform  of  [Vy(t,z=0)J  at  a  Fixed  x. 

Vy(co,z=L)  =  Fourier  transform  of  [Vy(t,z=L)j  at  a  Fixed  x. 
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V y(t,z)  is  a  voltage  defined  as  the  line  integral  of  Ey  from  the  microstrip  to  the 
ground  plane. 

Vv(t)  =  I  Ev(t,y)  dy,  at  fixed  x  and  z 

3 »  (5.10) 

where  the  integral  was  evaluated  by  using  the  4-th  order  Gaussian  quadrature.  Stable 

results  can  be  obtained  by  using  this  choice  of  voltage  definition  for  the  calculation  of 

the  frequency-domain  characteristics. 

The  effective  dielectric  constant  can  be  derived  from  the  calculated  phase 
constant,  (3,  by 

eeftfco)  =  P2(co)  /  o>2^0e0  (5.11) 


The  voltage-current  definition  is  used  for  the  characteristic  impedance 
calculation. 

Z(<o)  =  Vy(co,z=0)  /  Iz(co,z=0)  (5.12) 

where  the  current  is  found  by  a  line  integral  of  the  H-field  around  the  microstrip  line 
at  z=0. 

f  Hx(y=yb)dx+  f  H  (x=xj)dy  +  f  Hx(y=yjdx 

J°  Jy.  (5.13) 


5. 2. 3. 2.  Discontinities 

Figure  5.4  shows  a  general  two-port  discontinuity  with  a  uniform  microstrip 
line  connected  to  each  port,  where  two  reference  planes  and  two  observation  planes 
are  indicated.  In  order  to  obtain  the  frequency-dependent  scattering  parameters  of  the 
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Fig.  5.4  General  two-port  discontinuity  with  a  uniform  microstrip  line  connected 
to  each  port,  where  two  reference  planes  and  two  observation  planes  are 
indicated. 

discontinuity,  the  incident  wave,  Ey,inc>  and  the  reflected  wave,  Ey,ref>  should  be 
calculated  at  some  point  on  the  input  observation  plane  and  the  transmitted  wave, 
Ey,trs>  at  the  corresponding  point  on  the  output  observation  plane.  Then, 

S 1 1 (co)  =  [Vy,rcf(a),zi)/ Vyiinc(®,zi)]  /exp[-2yi(a>)zi]  (5.14) 

S2i(co)  =  [Vyitrs(co,z2)  VZoi(co)  /  Vy,inc(o>,z2)  VZ()2M  ]  /  exp[-yi(a>)zi- 
Y2(oo)z2)l  (5.15) 

where  yi(co),  Y2(a>)  are  the  propagation  constants  of  the  uniform  microstrip  lines 
connected  to  each  port  of  the  discontinuity.  Note  that  only  Si  j  exists  for  a  one-port 
structure,  such  as  an  open-end. 


5.3.  Results  and  Discussion 
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In  this  chapter,  four  structures  on  alumina  subtrate  (Ej-  =  9.6,  h  =  0.635  or  0.7 
mm)  shielded  by  a  rectangular  waveguide  (a  =  2,  b  =  10,  and  c  =  45  mm)  were 
studied.  These  are  a  uniform  microstrip  line,  a  step-in-width,  an  open-end  and  a  gap. 
In  the  calculation  of  the  field  variation  in  the  time  domain,  the  following  number  of 
discretization  points  were  used. 

number  of  discretized  points  in  x-direction  =  20  (Ax  =  0.095  mm) 
number  of  discretized  points  in  z-direction  =  100  (Az  =  0.445  mm) 
number  of  modes  =  19 

As  a  rule  of  thumb,  the  number  of  discretized  points  is  chosen  so  that  the 
discretization  length  is  approximately  one-tenth  the  length  of  the  smallest  wavelength 
to  be  considered. 


z 


(c)  (<n 

Fig.  5.5  Ey  configulation  beneath  a  uniform  microstrip  line  (a=2,  b=10,  h=0.635, 
W=0.635  mm)  (a)  t=0,  (b)  t=40,  (c)  t=80  and  (d)  t=  1 20  p-sec. 
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Figure  5.5  shows  the  pulse  propagation  along  the  uniform  microstrip  (h  = 
0.635,  W=0.635  mm)  in  the  time  domain.  In  order  to  check  the  convergence  of  the 
method  with  respect  to  the  number  of  modes,  the  effective  dielectric  constant  is 
calculated  using  a  different  number  of  modes.  As  shown  in  Fig.  5.6,  convergent 
results  can  be  found  by  using  more  than  seventeen  modes  for  the  tested  structure. 


Fig.  5.6  Convergence  of  the  method  with  respect  to  the  number  of  modes. 
A=19,  B=17,  C=15,  D=13,  E=1 1  and  F=curve  fit  formula  [36] 


Figure  5.7  shows  the  propagation  constant,  (3,  obtained  by  using  the  Fourier 
transform  of  the  time-domain  data  shown  in  Fig.  5.5.  The  results  agree  well  with 
those  from  the  curve-fit  formula  [36],  Figure  5.8  shows  the  calculated  transmission 


coefficient  which  is  supposed  to  be  one.  The  deviation  from  the  exact  value  is  less 
than  1  %. 


Fig.  5.7  Dispersion  relation  for  the  uniform  microstrip  line  (  a=2,  b=10,  h=. 635 
[mm],  e  1=1,  and  e 2=9.6)  obtained  by  TDML. 

0  TDML  ,  •  curve-fit  formula  [36]  (open  structure) 


Figure  5.9  shows  the  characteristic  impedance  (Zo)  of  the  structure  from  both 


the  TDML  calculation  and  the  closed  form  formula  given  in  reference  [36].  The 
difference  between  the  two  results  is  about  6  %  and  it  is  believed  to  be  caused  by  the 
side  wall  effect  which  is  known  to  reduce  the  characteristic  impedance  [37]. 


freq  [GHz] 


Fig.  5.9  Characteristic  impedance  for  the  uniform  microstrip  line  shown  in  Fig.  5.5. 
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Fig.  5.1 1  Scattering  parameters(Sl  1,  SI 2)  for  the  symmetric  step  discontinuity 


shown  in  Fig.  5.10. 


TDML, 


S 1 1  ref.  [26] 
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Figure  5.10  shows  the  scattering  of  the  pulse  at  a  microstrip  step  discontinuity 
in  the  time  domain.  Using  equations  (5.14)  and  (5.15),  the  frequency-domain 
characteristics  can  be  extracted  from  the  time-domain  data.  Figure  5.1 1  shows  the 
frequency-dependent  scattering  parameters  of  the  given  structure  with  published  data 
[26], 

Figures  5.12  and  5.13  show  the  scattering  of  a  pulse  at  a  microstrip  (h=0.7, 
W=0.7  mm)  open-end  and  a  gap  (s=0.35  mm)  in  the  time  domain,  respectively.  The 
frequency-dependent  scattering  parameters  derived  from  the  time-domain  data  are 
shown  in  Figs.  5.14  and  5.15.  Also  plotted  in  these  Figures  for  comparison  are  the 
FDTD  results  reported  in  [26].  The  results  agree  reasonablely  well  with  [26], 
although  they  do  show  some  undulation  due  to  discretization  errors.  It  is  believed 
that  this  error  can  be  reduced  by  using  a  non-uniform  discretization  scheme  without  a 
significant  increase  in  the  computation  time.  Finally,  this  method  can  be  easily 
extended  to  other  types  of  planar  transmission  structures  shown  in  Fig.  1.3. 

5.4.  Conclusion 

It  has  been  demonstrated  that  the  time-domain  method  of  lines  is  a  reliable  and 
efficient  method  which  is  capable  of  dealing  with  pulse  pmpagation  and  scattering  in 
discontinuity  problems  in  planar  trar  mission  lines.  This  is  an  extension  of  the  two- 
dimensional  results  described  in  the  previous  chapters.  It  has  been  found  that  the  DC¬ 
mode  should  be  included  in  the  modal  field  calculation  to  obtain  correct  results.  The 
field  can  be  described  accurately  by  the  proper  placement  of  the  discretized  Ey  field 
points  near  the  edge  of  the  microstrip  line.  Also,  the  frequency-domain  characteristics 
obtained  from  the  time-domain  data  have  been  made  more  stable  by  applying  the  QR 
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algorithm  in  the  eigenmode  calculation  and  by  using  the  integrated  field  quantities  in 
the  extraction  of  frequency-domain  parameters.  The  method  can  be  appiieu  to  other 
planar  structures. 


Fig. 5. 14  The  magnitude  and  phase  of  SI  1  for  the  microstrip  open-end 
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Chapter  6 
Conclusion 


6.1.  Achievements 

A  new  time-domain  method,  termed  the  Time-Domain  Method  of  Lines,  was 
proposed  for  the  analysis  of  microwave  and  the  millimeter-wave  planar  transmission 
lines  and  their  associated  discontinuities.  This  method  makes  use  of  the  Method  of 
Lines  technique  in  solving  the  time-dependent  Maxwell's  equations.  The  unique 
structural  feature  of  planar  circuits,  i.e.,  the  space  above  and  below  the  interfacing 
surface  is  uniform  and  homogeneous,  can  be  incorporated.  Therefore,  three- 
dimensional  structures  are  discretized  only  on  the  two-dimensional  substrate  surface 
and  an  analytical  solution  is  used  along  the  normal  direction.  Because  of  this  semi- 
analytical  property  of  the  method,  it  has  a  number  of  potential  advantages  over  many 
other  time-domain  methods.  First,  memory  and  computation  time  can  be  reduced 
since  analytical  processing  is  employed.  Second,  this  method  can  handle  an  open 
boundary  in  the  vertical  (y)  direction  easily  because  of  an  analytical  solution  used  in 
the  y -dire'' ’on. 

With  this  new  method,  the  accurate  cutoff  frequency  spectra  of  a  partially 
filled  rectangular  waveguide  and  a  finned  rectangular  waveguide  were  found  from  the 
impulse  response  of  the  structures.  Also,  the  method  was  extended  to  three- 
dimensional  problems  such  as  pulse  propagation  and  scattering  in  planar  transmission 
lines  with  discontinuities.  The  characteristic  impedance  and  the  propagation  constant 
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for  a  microstrip  line  and  the  scattering  parameters  of  its  discontinuities  (step-in-width, 
open-end  and  gap  discontinuities)  for  a  wide  range  of  frequencies  were  obtained  from 
the  time-domain  data. 

6.2.  Future  work 

In  future  research,  the  proposed  method  needs  to  be  applied  to  the  analysis  of 
other  discontinuities  in  many  planar  transmission  line  circuits.  The  radiation  effect  of 
a  discontinuity  can  be  analyzed  by  the  solution  of  an  open  boundary  problem  in  the  y- 
direction  by  using  the  proposed  method. 

The  method  itself  can  also  be  modified  to  improve  the  accuracy  of  the  results 
without  a  significant  increse  in  the  computation  time.  This  is  possible  through  the  use 
of  a  non-uniform  discretization  scheme  [33]. 

Since  most  of  the  computation  time  in  this  method  is  spent  on  the  calculation 
of  the  eigenfrequencies  and  the  corresponding  eigenmodes,  it  is  desirable  to  devise  a 
scheme  to  avoid  that  calculation.  The  two-region  approach  is  a  candidate  for  a  new 
scheme,  where  the  original  problem  is  divided  into  two  regions  (above  and  below  the 
substarte  surface)  and  solved  in  each  region  with  the  equivalent  magnetic  current 
source  on  the  interface  boundary  [38].  For  example.  Fig.  6.1  shows  the  equivalent 
two-region  problem,  where  the  original  inhomogeneous  problem  shown  in  Fig.  4.1  is 
divided  into  two  homogeneous  structures  with  equivalent  source  on  the  aperture 
region. 


M  =  n  X  E 


(6.1) 
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Fig.  6.1  An  equivalent  two-region  problem  of  the  probem  considered  in  Chapter  4. 

Also,  Fig.  6.2  shows  the  another  equivalent  structure  derived  from  Fig.  6.1 
applying  the  image  theorem  [39].  Notice  the  structure  6.2  is  a  simple  to  solve 
because  it  is  homegeneous  and  every  discretized  line  is  really  independent  each  other 
so  that  the  procedure  described  in  the  chapter  3  can  be  used.  In  this  approach,  the 
problem  is  to  find  the  accurate  equivalent  magnetic  current  on  the  aperture  at  each  time 
step.  It  can  be  obtained  by  applying  the  finite  difference  scheme  around  the  aperture 
region  in  Fig.  6.1. 


Fig.  6.2  An  equivalent  structure  of  Fig.  6.1  using  the  image  theorem. 


_N+1..  ..  „N..  ..  _N-L.  .. 

Ez  (i,j)  —  2Ez(i,j)  ■  Ez  (i,j) 


+  ^-jE^(i+l  j)  -  2E^(i  j)  +  Ej(i-l,j)) 


c  At  In 


(Ez(i,j+l)-2Ez(io)  +  E^i,j-l) 


where  superscript  N  represents  the  N-th  time  step. 


The  difficulties  encountered  in  this  approach  were  in  the  representation  of  the 
pulse  in  terms  of  continuous  basis  functions  on  each  line  accurately  and  in  the 
reproduction  of  original  wave  in  terms  of  sampled  pulses. 


Appendix  1 

Determination  of  the  Eigenvalues  and  Eigenvectors  of 
Second-Order  Difference  Matrix 


The  diagonalization  matrix  can  be  obtained  from  the  eigenvectors  belonging  to 
the  eigenvalues  of  the  real  symmetric  matrix  [D^j. 


([Dxx]-Jik[I])[X]k  =  0 


(All) 


where  [Dxx]  is  given  by  equation  (2.14)  and  the  subscript  k  means  the  k-th 
eigenvalue  and  eigenvector. 

Since  [Dxx]  is  a  tridiagonal  matrix,  equation  (A  1.1)  actually  represents  N  second- 
order  difference  equations 


00 


(k)  v  (k) 


-x;x2->.k)x','-x;;;  =  o 

>1  -  1  , 


N 


(A  1.2) 


where  the  subscript  i  means  the  i-th  element  of  the  k-th  eigenvector. 
Using  the  solution  of  the  form, 

X  :k>=  A ,  e’,<Pt+  B  .  e 


(A1.3) 


the  characteristic  equation  can  be  obtained  from  (A  1.2),  which  has  the  characteristic 
root 

2 

X  k  =  2  (  1  -  cos  <p  k ) 


(A  1.4a) 


or 


A. .  =  4  sin  <p  k  /  2 


7  2 


(A1.4b) 
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where  tp^  depends  on  the  boundary  conditions  and  are  determined  by  using  the 
boundary  conditions  shown  in  Table  2.1. 


(i)  Neumann  -  Dirichlet  Condition 

Assuming  the  number  of  the  discretized  lines  is  equal  to  N,  the  first  and  the 
last  equation  become 


x0(k')  -  X,,k)  =0 

XN+,(k)=o 


or 


1  -  e  ‘J<Pk 

e  -i(N+  1  )  <pk 


(A  1.5a) 
(A  1.5b) 


(A  1.5c) 


In  order  to  have  nontrivial  solutions, 

w  _  k  -  1/2  ,  ,  .  XT 

Mj.  i  /  ?  ^  —  1,2,...,N 

N+1/2  (A  1.6) 

Using  the  first  row  of  equation  (A1.5)  and  the  equation  (A1.3), 


(  1  -  eJ<Pk)Ak=-(  1  -  e'J<Pk)Bk 
e J<Pk/  2  Ak  =  e  J(Pk/  2  Bk 
x!k)=  Akcos(i  -  l/2)<pk 


(A1.7a) 
(A  1.7b) 
(A  1.8) 


After  normalization,  the  formula  for  the  elements  of  the  transformation  matrix  can  be 


obtained. 
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ND.ik  ~  y 


N+l/2 


cos 


(i-  l/2)(k-  1/2)* 
N+l/2 


i,k=  1,2 . N 


(A  1.9) 


(ii)  Dirichlet  -  Neumann  Condition 

Assuming  the  number  of  the  discretized  lines  is  equal  to  N,  the  first  and  the 
last  equation  becomes 


XqOO  =  0 

-  xN(k)  +  XN+100  =  0 


or 


cjN«.(eW._i)  e  .,N<p.(e  -rn_  I} 


(A  1.10a) 
(Al.lOb) 


(Al.lOc) 


In  order  to  have  nontrivial  solutions, 

k  -  1  /  2  ,  .  _  vT 

<Pk_N+  1  12  K’  k"  1’2’-'-N 
Using  the  first  row  of  equations  (Al.lOc)  and  (A  1.3), 


(A  1.11) 


Ak  =  *Bk  (A1.12a) 

X  i  -  Ak  sin  i  (pk  (A  1.12b) 

After  normalization,  the  formula  for  the  elements  of  the  transformation  matrix  can  be 
obtained. 


T 


ND.  ik  ~  'y 


N  +  l/2 


sin 


i  (  k  -  1  /  I. )  k 
N  +  l/2 


i,  k  =  1 , 2, ...,  N 


(A  1.13) 


(iii)  Dirichlet  -  Dirichlet  Condition 
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Assuming  the  number  of  discretized  lines  is  equal  to  N,  the  first  and  the  last 
equation  become 


X0(k)  =  0 
XN*i(k>  =  0 


or 


i(N+  1  )<pk  -j(N+  1  )<pk 

V  V 


Ak 


=  0 


In  order  to  have  nontrivial  solutions, 

^k  =  jq  +  i ’  k=U2,...,N 
Using  the  first  row  of  equations  (Al.l-c)  and  (A1.3), 


(A1.14a) 

(A1.14b) 


(A1.14c) 


( A 1 . 1 5) 


v  (k)  .  .  . 

X  j  =  Ak  sin  l  9k 


( A 1 . 1 6) 
(A  1.17) 


After  normalization,  the  formula  for  the  elements  of  the  transformation  matrix  can  be 
obtained. 


1 


DD.  ik  -  Y 


N  +  1 


sir 


i  k  7i 

nTT 


i,k  =  1,2 . N 


(A  1.18) 
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(iv)  Neumann  -  Neumann  Condition 

Assuming  the  number  of  the  discretized  lines  is  equal  to  N+I  to  be  consistant 
with  the  D-D  case,  the  first  and  the  last  equation  become 


XqOO  -  Xj(k)=  0 

>wk)  -  o 

d-eJ®')Ak 
<l-e'*‘)Bk 

In  order  to  have  nontrivial  solutions. 


1  1 

j  (  N+  1  )  <pk  -j  ( N  +  1  )  <p. 
e  e 


(A1.19a) 

(A1.19b) 


(A1.19c) 


k  7t 

(Pk  =  NTT’  k  =  0’ 

Using  the  first  row  of  equations  (A  1.1 9c)  and  (A  1.20), 

(  1  -  eJ<Pk)  Ak  =  -  (  1  -  e  J<Pk)  Bk 
eJtPk/2Ak  =  e-J<Pk/2Bk 
X  -k)=  Ak  cos  ( i  -  1  /  2  )  <pk 


(A  1.20) 

(A1.21a) 
(A  1.2  lb) 
(A  1.22) 


After  nomialization,  the  formula  for  the  elements  of  the  transformation  matrix  can  be 
obtained. 


T 


NN 


w 


1 


N  +  1 


T 


NN,  i  k 


W 


N  +  1 


.os 


(i-  1  /  2  )  k  re 
N  +  1 


i  =  1,2 . N+ 1 ,  k  =  1,2 . N 


(A1.23) 


Appendix  2 

Calculation  of  the  Quasi-Diagonal  Matrices  [dx]  of  the  First- 
Order  Difference  Operator  [Dx] 

The  matrix  [dx]  is  defined  as  a  product  of  the  first-order  difference  matrix 
[Dx]  with  the  corresponding  transformation  matrix  [Tx]  from  the  right  and  with  the 
transposed  matrix  [Txd]t  of  the  dual  boundary  value  problem  from  the  left. 

[dx]  =  [Txd]t  [Dxl  [Tx]  (A2.1) 

where  the  superscript  'd'  denotes  the  dual  boundary  condition.  Before  [dx]  is 
determined,  one  result  may  be  anticipated.  The  eigenvalue  matrix  [X.2]  results  from 
the  product 

[dxl1  [dxl  =  Pxl1  [Dxlr  [Txdl  rTxdll  [Dx]  [Tx]  =  [Tx]t  [Dx]t  [Dx]  [Tx]  =  [\2] 

(A2.2) 

w'here  the  equations  (2.16)  and  (2.14)  are  used.  Therefore,  the  relation  between  [dx] 
and  [A-]  is  anticipated. 

(i  )  Boundary  Conditions  DD  and  NN 

For  the  boundary  condition  DD.  the  first-order  difference  matrix  [Dxdd]  is 
given  in  Table  2.2.  The  coiresponding  [dxDD]  is 

[dxDDl  =  \JxW]'  [DXD^1  [TXDD] 


7  7 


(A2.3) 
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If  we  multiply  [TxNN]t  by  [DXDD],  we  find  that  the  first  row  of  the  result  vanishes 
because  all  the  elements  of  the  first  column  vector  in  [TXNN]  are  the  same.  The 
remaining  part  of  the  resulting  matrix  can  be  written  as  a  product  of  and 

[Txdd],  where  [aDl>]  is  the  diagonal  matrix  foimed  by  the  positive  square  roots  of 
[>vDd  23-  Therefore, 


T 


NN 


.it 


D,1  = 


0 


•DO 


DD 


(A2.4) 


Since  [TXDD]  is  a  symmetric  matrix,  [TXDD]  =  [Tx130]1.  Consequently  [dxDD]  is 
given  by 


0 

^  DD 


Thus  [dxDD]  is  a  quasidiagonal  matrix. 

Using  the  equation  (2.5),  it  can  be  shown  easily  that 


(A2.5) 


[dx™]  =  -  [dxDD] 


(A2.6) 


for  the  boundary  condition  NN. 


(ii)  Boundary  Condition  DN  and  ND 

For  the  boundary  condition  DN,  the  first-order  difference  matrix  is  given  in 
Table  2.2.  In  tliat  case,  the  corresponding  fdxDN'j  is  given  by 


idxDN1  =  (V^)' fDXDNHTxDN] 


(A2.7) 
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Using  the  procedure  described  previously,  we  can  find  that 

[dxDN]  =  [XDN]  =  [dxND]  -  (A2.8) 

where  the  k-th  diagonal  element  of  [^dn1»  ^-k>  given  by  the  positive  root  of  X^-. 


Appendix  3 

Calculation  of  Modal  fields  in  three-dimensional  problems 

(I)  Time-Harmonic  Case 

In  three-dimensional  problems  shown  in  Fig.  5.1,  all  the  six  Field  components 
should  be  considered  to  satisfy  the  boundary  conditions.  Usually,  two  scaL' 
potentials,  lFe(x,y,z,t)  and  vF^1(x,y,z,t),  are  introduced  to  simplify  the  formulation 
instead  of  considering  all  six  Field  components.  The  Fields  can  be  derived  from  the 
solution  of  two  scalar  potentials.  The  equations  for  the  potentials  are  given  by 

32  /  5x2  +  52  vpe  /  3y2  +  32  ye  /  3Z2  -  32  vpe  /  3t2  =  0  (A3. 1  a) 

32  q/h  /  ax2  +  52  q/h  /  3y2  +  32  vj/h  /  3Z2  -  32  vph  /  at2  =  0  (A3.1  b) 

Using  the  separation  of  variables  technique,  'Ffx^z.t)  =  y(x,y,z)  T(t),  the 
time-dependent  solution  is  of  the  form 

T(t)  =  A  cos(cot)  +  B  sin(cot).  (A3. 2) 

The  space-dependent  equations  are  given  by 

52  /  5x2  +  32  ye  /  3y2  +  52  /  dz2  +  co2(ie(v)  y^  =  0  (A3. 3a) 

52  yh  /  5x2  +  32  yh  /  3y2  +  32  xj/h  /  3z2  -r  co2fie(y)  y^  =  0  (A3. 3b) 

with  the  boundary  conditions. 
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V^x.y,*)  =  0,  at  y  =  0,  b 
9yh(x,y,z)  /  9y  =  0,  at  y  =  0,  b 


(A3.4a) 
(A3. 4b) 


At  y  =  h. 


2  e  h  2  e  h 

3  Vi  9  Vi  a  vn  3  Vn 

Exi-Exn* - * - +  =  0 


jOEjEgSx  dz  dy  jcjJEgBx  dz  dy 


(A3.4c) 


1  92  2  e  1  92  2  e 

Ezi-  Ezn= - (— +  erco  ^q)  y, - (— +co 

JClEjEo  dz  jox0  dz 


(A3.4d) 


<■)  h  e  e 


=  _j_ivi  aVi  i  a  vn  ayn_ 

xI  xI1  jco|i0  ax  9z  +  9y  j(0(Io  9x  3z  9y 


(A3.4e) 


1  92  2  h  1  92  2  h 

Hzi-HzII  = - (— +  0)  M-oEjEp)  Vj - (— -  +  co  ^o^o)  V,i=  J, 

Mi0  dz  jc4i0  dz 


(A3.4f) 


where  subscripts  I,  II  represent  the  regions  shown  in  Fig.  5.1. 


The  first  step  of  analysis  is  the  discretization  of  the  equation  according  to  the 
rule  explained  in  Chapter  2,  which  gives 


e  ND  e  e  ^  ND'  _  I 

C*  11/  Dxx  V  l)/  D.x  2  e  r, 

— ^2  *  J2  ~  +  rJLj>!-±+  co  My)  [v  J  =  [oj 

ay  Ax  Az 


(A3. 5a) 
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(  jPr.TDN][v]  Jv]lP»x  J-4-  CoWy)[v]  =  [°] 


(A3.5b) 


with  the  boundary  conditions 
[yC]  =  0,  at  y  =  0,  b 
9[Yhl  /  dy  =  0,  at  y  =  0,  b 


(A  3. 6a) 
(A3. 6b) 


At  y  =  h. 


7T  iTii  3y  j^o  AxAz  uy 


(A3. 6c) 


1  Vi 

_ ! (L—Ud — - 

v  2 
jcce^o  Az 


]4  -  jl  dyj& + 444  ■  W 

W  a,2  1 


+  erco  ^^dVij) 


jcoe0  Az 


(A3.6d) 


pnjaPj 

jco^o  Ax  Az 


i'I4K1  41 

+  3y  '»  A*a*  y 


(A3.6e) 


hlr~DNf 


h|fnDN' 
1  Vl  j  1‘^lZ 


j_  <l±4 + e  444  - — + “  v#w = p  J 

'  .  2  r  ^  j£4i0  Az 

Jto^o  Az 

(A3. 60 


Since  the  matnces  ID,,™].  ID*,""].  Pzz^l.  and  [D*z«*l  a- 
symmetric  matnces,  they  can  be  dtagonaHzed  by  appropriate  mansfomtarion  mamces 


(see  Table  2.3).  The  potentials  can  also  be  transformed  into  the  "transformed 
potentials". 


[U]  =  [Tx^]1  m  [T^] 

[V]  =  [TxDN]t  [vh]  [T^Nj 

Then,  equations  (A3.6a,b)  are  transformed  into 


aiu]  + ' 

'  ND 

dxx 

[m+m 

’  ND 

dZ2 

1_ 4.  /-.x  1  ip/v\  Ft il  —  f ol 

By2 

Ax' 

'  DN 
dxx 

) 

Im  ^  m( 

2 

Az 

'  DN 
dzz 

co  p.tty;  [uj  — [uj 

t 

_ /•.a  r\/i-  f nl 

ay2 

2  2  ^  ;  L  J  —  L  J 

Ax  Az 

with  the  transformed  boundary  conditions, 
[U]  =  0,  at  y  =  0,  b 
3[V]  /  9y  =  0,  at  y  =  0,  b 


At  y  =  h, 


(A3.7a) 

(A3.7b) 


(A3. 8a) 

(A3. 8b) 

(A3. 9a) 
(A3.9b) 


(A3.9c) 

[o] 
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1 


dr]'fv.l[df]  ^[U.l  1  [dfffVll] 

oy  " ' 


,ND' 


jcojj.0  Ax  Az 


Ax  Az 


U 

57 


^•[ji 


(A3.9e) 


r^iTjDNl  „  r-.7  iLDN] 

1  /[^l][^zzj  2  1  JVIl][dzz  2  r, ,  Cf-] 

—  (L  J  2  J-  +  co  ^(JVj)  -- - (L  -I1  2  +(0  HdE(|Vn])=[jxj 


jcoM-o  Az‘  '  '  jcop-o  Az‘ 

(A3.90 

Now  equations  (A3.8a,b)  are  uncoupled  so  that  the  general  solution  of  the  ik- 
th  line  which  satisfies  the  boundary  condition  (A3.8a,b)  can  be  obtained 


UUk=  Ai,lksinh  K,  iky 
Uii,ik=AIItiksinhKIUk(b-y) 

vi. 1k=Bi,ikcoshTii,iky 

vii, ik=Bn,ikCoshqIIiik(b-y) 


(A3. 10a) 
(A3. 10b) 
(A3. 10c) 
(A3.10d) 


where 

Ki.ik2  =  -  [dxx.i^  /  Ax2  +d2ZjkND  /  Az2  +  coV^l  (A3. 1 1  a) 

Kn.ik2  =  *  [dxx,^  /  Ax2  +dZZrkND  /  Az2  +  co2p£o]  (A3. 1 1  b) 

rji,ik2  =  -  [dXx,iDN  /  Ax2  +dzz>kDN  /  Az2  +  co2q£r€o]  (A3. 1 1  c) 

Tln.ik2  =  -  [dxx,iDN  /  Ax2  +dZZ;kDN /  Az2  +  co2q£o]  (A3.1  Id) 


By  means  of  these  solutions,  the  derivatives  of  U 
h  can  be  represented  by  the  values  of  U  and  V  at  y  =  h. 


dU 


I.ik 


‘I.ik 


dy 

dUn.ik 


tanh  Kj  ik  h 


Ui.ik 


‘II, ik 


dy 


tanhKnilkd 


,  at  y  =  h 

U,i,lk 


,  at  y  =  h 


and  V  with  respect  to  y  at  y 


(A3. 12a) 

(A3. 12a) 
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dV  i.ik 
dy 

dV  II. ik 
dy 


=  7li..ktanhTli.ikhVi.ik  t 

,  at 

=  -Tll.iktanhTlukdVII..k 


(A3. 12a) 
-  (A3. 12a) 


Now  the  equations  (A3.9c~f)  have  six  unknowns  Uiiik>  Un^k,  Vjik,  Vn,ik> 
Jx  ik  and  Jz jk  at  y  =  h.  Therefore,  the  continuity  equations  in  (A3.9c~f)  can  be  used 
to  solve  for  the  potentials,  U^ik,  Un,ik,  Vi jk  and  Vq^k  at  y  =  h  in  terms  of  the  current 
densities,  Jx,ik  and  Jz-±. 


U,' 

1 

c22  h  l  ‘  cl2d2 

C1 2 

Jx 

Vi. 

ik  c  1 1  c22 “ c 1 2  c2 1 

cub2- c2i  bi 

‘  C1 1 

(A3. 13a) 


where 
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jcoe^Q  Ax  Az 
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Az 
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c i  2  =  - 11 1,ik  tanh  q Itlk  h  -  T1  Iiilk  tanh  q  JI  ik  d 
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Az 


+  w  Wo 
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+  C0  n  of0 


Az 


(A3. 13b) 


(A3. 13c) 
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/  d1®  2  \ 


c21  =  KIiikCOth  KTj  ikh  +  KIItikCOth  KII>ikd 


Az 
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d77  k  ^ 
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Az 
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,DN  ,DN 
1  ^x,i  ^z,k 


jcoji0  Ax  Az 
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Az 
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dzz,k  2 
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Az 


+  “¥o 


b.  =-jco^0TiIIiktanhriIIikd 


d™..  2 


lzz,k 


Az 


+  CO 


b9  =  - 


,DN  ,DN 
dx.i  dz.k 
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Ax  Az  2 

uzz,k 


Az 


/  (A3. 13d) 


(A3.l3e) 


(A3. 13f) 


(A3.13g) 


Then,  the  field  values  at  y  =  h  can  be  found  by  using  the  following  relations  evaluated 
at  y  =  h. 


,ND  ,ND 
1  dx,.dz,k 


Ey,.i/y)  = 
Ez.!k(y)  = 


jcoe(y) 

Ax  Az 

1 

,ND 

dz,k  a  . 

jcoe(y) 

Az  ^ 

1 

,ND 

Ulk(y)-Tllk^nhriikh  Vlk(y) 


,DN 

dx.. 

Ax 


(A3. 14a) 

(A3. 14b) 


jcoe(y) 


Az 


(A3. 14c) 
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H*.ik(y)  = 


1C; 


,DN  ,DN 

iV  1  0.  ;  Cl 7  ^ 

Uik(y)+— - Sil-£jiVlt(y) 


tanh  Kikh 

X,1 


jo^M-o  Ax  Az 


Hy.ik(y)  =  -  —  ulk(y)  + 
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1  dz,k  0 
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jtop.0  Az 
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Wo 


Az 


(A3.14d) 

(A3.14e) 

(A3. 140 


From  this  point  on,  we  change  the  two-dimensional  arrays  for  fields  and 
currents  into  one-dimensional  arrays  for  convinience  of  formulation  according  to  the 
rule 


n  =  i  +  (k  -  1)  Nx  ,  1<  i  <  Nx,  1<  k  <  Nz 


(A3. 15) 


Then,  tne  two-dimensional  Nx  by  Nz  matrices  become  NXNZ  column  matrices  and 
equations  (A3.14a,b)  can  be  expressed  in  the  following  form: 


(A3. 16) 


where  [Z]'s  are  diagonal  matrices. 

Before  the  final  boundary  condition  is  applied,  a  useful  matrix  operation,  the 
Kronecker  product,  is  introduced 


[A]  ®  [B]  = 


31  lP] 

am{B] 


aUBJ 

amr[B] 


(A3. 17) 
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Two  important  properties  of  the  Kronecker  product  are 

( [A]  ®  [B]  )  ( [C]  ®  [D] )  =  ( [A]  [C] )  ®  ( [B]  [D] )  -  (A3. 18a) 

( [A]  ®  [B]  )t  =. [A]1  ®  [B]t  (A3. 18b) 


Now,  the  final  boundary  condition  is  that  the  tangential  electric  fields  should 
be  zero  on  the  metalization  part  of  the  substrate  surface.  In  order  to  apply  this 
condition,  we  have  to  go  back  to  the  original  domain  using  the  inverse  transformation 
relation 
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ND 
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IP 
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— >  [j  x] = ( p 


,CN 


® 
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Tx 

NDl 


)[eJ 


[jj=([T^]®[Tr])frj 


(A3. 19a) 
(A3. 19b) 
(A3. 19c) 
(A3.i9d) 


where  the  subscript  2D  means  the  two-dimensional  matrix  considered  up  to  equation 
(A3. 14).  The  equation  (3.15)  can  be  inverse  transformed  into 


rE  n 

'  N' 

N 

N 

i- .j 
_ 1 

A 

A 

red 

p-xzjp-xx 

red 

A 

(A3. 20) 


where  red'  represents  the  reduced  matrix  made  by  deleting  the  rows  and  columns 
corresponding  to  the  non-metalization  lines  from  the  full  matrix.  Equation  (A3. 20) 
will  have  nontrivial  solutions  at  the  resonant  frequencies  of  the  structure  only.  The 
resonant  frequencies  are  found  from  the  determinant  equation 


8  9 


det 


'Mz»f 


(A3.21) 


A:  resonance,  the  modal  current  distribution  can  be  found  as  an  eigenvector  in 
(A3.20)  by  using  the  QR  algorithm  for  the  stability  of  the  solution. 

One  inportant  thing  to  be  mentioned  in  the  calculation  of  the  roots  is  the 
elimination  of  poles  present  in  the  matrix  element  in  order  to  obtain  the  correct  roots. 
One  way  to  remove  the  poles  in  the  determinant  calculation  is  by  the  pole 


multiplication  method,  in  which  all  the  elements  are  multiplied  by  the  common 
di.iominator.  In  this  case,  the  pole  extraction  is  performed  in  equation  (A3. 16),  not  in 
equation  (A3.21).  This  is  because  the  poles  in  equation  (A3. 16)  can  be  easily  found. 
The  process  is  as  follows: 


where 

N 

Di=  n  Dj 

i=l,i*i 


(A3. 22b) 


After  the  values  of  the  potentials,  Ui,ik,  Un,ik>  Vitik,  Vn,ik  at  y  =  h  are  known 
from  the  modal  current  density  distribution  obtained  in  the  previous  step,  the 
coefficients  Aj  ^,  Bj  ^,  Ajjjk  and  B^ik  in  equations  (A3.10a~d)  can  be  determineo. 

Therefore,  the  transformed  modal  fields  of  the  given  structure  can  be  obtained  by 
equations  (A3.14a~f). 
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(II)  Static  Case 

In  this  case,  the  governing  equation  is  the  Laplace's  equation  with  proper 
boundary  conditions: 


2 

V  V-0 
and  at  v  =  h, 


av,  avn 

E^'E^~"57+_5r"0 


E,n  =  - 


av,  avn 


'zi'  czn-  --V— +  -a —  =  0 


avj  avn 

erEyi-Eyn--er-^r+-^=-P 


(A3. 23) 

(A3. 24a) 

(A3.24b) 

(A3. 24c) 


where  the  subscripts  I  and  II  represent  the  regions  shown  in  Fig.  5.1. 

Using  the  discretization  ru'e  described  in  Chapter  2,  equations  (A3. 23)  and 
(A3.24a~c)  can  be  discretized  as  follows  : 


ay4 

at  y  =  h 
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ND 
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Df][V„]=0 


D?j-0 


(A3. 25) 


(A3. 26a) 


(A3. 26b) 


(A3. 26c) 
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fw-[p] 

dy 


where  [  ]  represents  a  matrix  whose  ij-th  element  represents  the  potential  at  the  ij-th 
line. 

Since  the  second-order  difference  matries  [DXXND]  and  [Dzzdn]  are  real 
symmetinc  matrices,  they  can  be  diagonalized  by  using  the  transformation  matries  (see 
Table  2.3) 

[Dxx^]  [TxND]  =  diag  [d^]  (A3.27a) 

[TzDN]t  [DzzDN]  [TzDN]  _  diag  [d^DN]  (A3. 27b) 


The  potential  [  V]  is  also  transformed  into  a  transformed  potential  [U], 


[U]  =  [TxND]t  [Y]  [TzDN] 


(A3.28) 


Then,  equations  (A3. 23)  and  (A3.24a~c)  become 


ay 

at  y  =  h 


^-[u]+^[d^][u]+-Ltu][d^' 

Ax  Az 
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(A3. 29) 
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(A3. 30a) 
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(A3. 30b) 


•erJ[U,]+i-[Un]  =  . 


(A3. 30c) 
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The  solution  of  equation  (A3. 29)  which  satisfies  the  top  and  bottom  boundary 
condition  can  be  written  as 


UI,ik=  AI,ik 


un,ik=  An,ik 


sinh  Kiky 

K 

lk  ,  in  region  I 
sinh  Kik(  b  -  y  ) 


kik  ,  in  region  II 


(A3. 31a) 


(A3. 31b) 


where 


2  DN 

^  _  dxx,i  d2zk 

Kik  2  '  2 

Ax  Az 


(A3. 31c) 


By  means  of  these  solutions,  the  derivative  of  the  potentials  with  respect  to  y 
at  y  =  h  can  be  represented  in  terms  of  the  potential  at  y  =  h. 


auIik 

~a7~y=h=  K.kcoth  Kikh  ui,ik(y=h) 


3Utt 


y=h 


=  -  Kikcoth  KikdUIliik(y=h) 


(A3. 32a) 


(A3. 32b) 


Now,  the  boundary  conditions,  (A3.30a~c),  can  be  solved  for  [Uj  jjJ  and 
[Un,ikl  at  y  =  h  in  terms  of  the  charge  density  [p]. 


Ui,ik(y=h)  =  UIIJk(y=h)  = 


erKikCOthlCikh  +  KikCOthKikd 


(A3. 33) 


Then,  the  field  values  at  y  =  h  can  be  found  by  using  the  following  relarions  evaluated 
at  y  =  h. 
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dNE> 

Ex.lk(y)  =  --^Ulk(y) 

Ax  (A3. 34a) 

Ez,ik(y)  =  -  — uik(y) 

Az  (A3. 34b) 


From  this  point  on,  the  two-dimensional  field  arrays  are  changed  into  the 
corresponding  one-dimensional  vectors  for  convenience  of  formulation  according  to 
the  rule  described  in  equation  (A3. 15).  Then,  the  two-dimensional  Nx  by  Nz  matrix 

becomes  an  Nx*Nz  column  vector  and  the  relations  (A3. 34a  and  b)  can  be  expressed 
in  the  matrix  form: 


(A3. 35) 


where  the  [Gxp]  and  [G2p]  are  diagonal  matrices. 

The  final  boundary  condition  is  that  the  tangential  electric  field  should  be  zero 
on  the  metal  strip.  In  order  to  apply  this  condition,  we  have  to  go  back  to  the  original 
domain  using  the  inverse  transformation  relation 


fE*H 


[EzH 


f  ON' 

[Tf])fEj 

tND 

1  z 

® 

1 

Tf] )  [EJ 

(A3. 36a) 
(A3. 36b) 


® 


(A3. 36c) 


where  <8>  is  the  Kronecker  product  defined  by  (A3. 17). 
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Then, 

[p]-  =  [0] 

red  (A3. 37) 

where  the  subscript  'red'  denotes  the  reduced  matrix  obtained  by  deleting  the  rows 
and  columns  corresponding  to  the  non-metallic  discretization  points.  The  static 
charge  distribution  can  be  calculated  from  equation  (A3. 37;  by  using  the  QR 
algorithm.  After  the  values  of  the  potentials,  Uj  ^  and  Ujj^k  at  y  =  h,  are  found 
from  the  static  charge  distribution,  the  coefficients  Aj  ^  and  Ajj  ^  in  (A3.31a  and  b) 
can  be  obtained  by  evaluating  the  potential  values  at  y  -  h,  from  which  the  electric 
field  distributions  can  be  calculated  by  equations  (A3. 34a  and  b). 
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